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ABSTRACT 

We investigate pathways for the formation of icy super-Earth mass planets 
orbiting at 125-250 AU around a 1 M 0 star. An extensive suite of coagulation 
calculations demonstrates that swarms of 1 cm to 10 m planetesimals can form 
super-Earth mass planets on time scales of 1-3 Gyr. Collisional damping of 
10 2 —10 2 cm particles during oligarchic growth is a highlight of these simulations. 

In some situations, damping initiates a second runaway growth phase where 
1000-3000 km protoplanets grow to super-Earth sizes. Our results establish the 
initial conditions and physical processes required for in situ formation of super- 
Earth planets at large distances from the host star. For nearby dusty disks in 
HD 107146, HD 202628, and HD 207129, ongoing super-Earth formation at 80- 
150 AU could produce gaps and other structures in the debris. In the solar 
system, forming a putative planet X at a < 300 AU (a > 1000 AU) requires a 
modest (very massive) protosolar nebula. 

Subject headings: planetary systems - planets and satellites: formation - solar 
system: formation 


1. INTRODUCTION 

For nearly a century, planet X has haunted the outer Solar System. Originally envisioned 
as a way to explain irregularities in the orbits of Uranus and Neptune (e.g., Lowell 1915; 
Bower 1931; Brouwer 1955; Rawlins 1970, 1972; Harrington 1988), planet X - along with the 
variants Nemesis and Tyche - has since been invoked to explain the orbital properties of some 
comets and trans-Neptunian objects (e.g., Matese et al. 1999; Gomes et al. 2006; Lykawka 
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& Mukai 2008) and apparent periodicities in mass extinction events on Earth (Davis et al. 
1984; Whitmire & Jackson 1984). Generating these phenomena requires a massive planet X, 
~ a few Earth masses to tens of Jupiter masses, with a semimajor axis a ~ 200-20,000 AU. 

Observational and theoretical constraints severely limit possibilities for planet X (e.g., 
Hogg et al. 1991; Zakamska & Tremaine 2005; Iorio 2009; Fernandez 2011). Along with 
improved mass measurements for Uranus and Neptune from Voyager, better analyses of 
positional data eliminate the need for a massive planet to explain any ‘irregularities’ in 
the orbits of Uranus and Neptune (e.g., Standish 1993). By excluding the possibility of a 
brown dwarf anywhere in the Solar System, a Jupiter-mass planet inside ~ 25,000 AU, and 
a Saturn-mass planet inside ~ 10,000 AU, observations from the IRAS, 2MASS, and WISE 
all-sky infrared surveys specifically rule out Nemesis and some other variants of planet X 
(e.g., Luhman 2014). 

The recent discovery of 2012 VP 113 (Trujillo & Sheppard 2014) in the ‘inner Oort cloud’ 
has renewed interest in planet X. With orbital parameters similar to those of Sedna 1 , 2012 
VP 113 joins a group of a dozen relatively large (R ~ 200-1000 km) distant objects ( q p > 
30 AU, e > 0.7) with arguments of perihelion 00 ~ 0°. The observed distribution of 00 for 
these objects is statistically unlikely (see also de la Fuente Marcos & de la Fuente Marcos 
2014). Although dynamical interactions between the young Sun and a passing star can 
produce objects on highly eccentric orbits like Sedna and 2012 VP 113 (e.g., Morbidclli & 
Levison 2004; Kenyon & Bromley 2004c), the known planets randomize u on 0.1-1 Gyr time 
scales (Gomes et al. 2006; Lykawka & Mukai 2008; Trujillo & Sheppard 2014). Super-Earth 
mass planets at a ~ 200-300 AU can maintain the observed distribution of 00 for the age of 
the Solar System (Trujillo & Sheppard 2014; de la Fuente Marcos & de la Fuente Marcos 
2014; Iorio 2014). 

In addition to explaining the orbital dynamics of objects beyond 40-50 AU, identifying 
one or more Earth-mass planets with a ~ 200 AU has interesting implications for the archi¬ 
tecture of the Solar System beyond Neptune. Despite a lack of massive planets with a > 
35 AU, there is a diverse set of trans-Neptunian objects with radii R < 1000-2000 km on 
orbits with a broad range of e and i (e.g., Gladman et al. 2008; Petit et al. 2011). Aside 
from objects like Sedna and 2012 VP 113 , the distribution of semimajor axes for these objects 
suggests a clear ‘edge’ at roughly 48 AU. Confirming a massive planet at ~ 200 AU might 
imply a very large reservoir of material beyond this edge. 

Establishing the existence of a massive planet beyond 100 AU has profound consequences 
for the origin of the Solar System. In current theory, planets grow out of material orbiting 
within a circumsolar disk of gas and dust. As the gaseous disk dissipates, these planets (i) 


1 Sedna: semimajor axis a = 542 AU, eccentricity e = 0.86, and perihelion distance q p = a( 1 — e) = 
76.23 AU; 2012 VP n3 : a = 266 AU, e = 0.698 AU, and q p = 80 AU 
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clear out material along their orbits and (ii) migrate radially to a more stable orbit (e.g., 
Youdin & Kenyon 2013; Morbidclli 2013). If an extra super-Earth forms with the known 
gas giants at a ~ 10-20 AU before the gaseous disk disappears, this planet could migrate 
radially to much larger semimajor axes as in models for Fomalhaut b (Crida et ah 2009). 

Planet-planet scattering is a promising alternative to radial migration (Bromley & 
Kenyon 2014). In this picture, several massive planets form out of disk material at 5- 
10 AU. Close orbital encounters among these planets then scatter lower mass planets to 
large semimajor axes (e.g., Rasio & Ford 1996; Levison et al. 1998; Juric & Tremaine 2008; 
Bromley & Kenyon 2011a). If the gaseous component of the disk is massive and extends to 
125-250 AU, gas drag and dynamical friction can circularize the orbits of super-Earth-mass 
planets on time scales comparable to the disk lifetime. Because lower mass planets are left 
on eccentric orbits, measuring the orbital elements of any planet with a > 100 AU constrains 
the mass and outer radius of the solar nebula (Bromley & Kenyon 2014). 

Establishing the orbital eccentricity of a super-Earth at 200 AU provides a way to 
distinguish between these two possibilities. Migrating planets have nearly circular orbits 
(e.g., Ward 1997; Masset & Papaloizou 2003; Ida & Lin 2008; Crida et al. 2009). Although 
gas drag can produce nearly circular orbits among scattered massive planets, lower mass 
planets have more eccentric orbits (Bromley & Kenyon 2014). 

To provide a distinct set of predictions for massive planets at 200 AU, we consider 
models where one or more super-Earths form directly out of solid material well beyond the 
outer edge of the Kuiper belt (see also Stern 2005; Kenyon & Bromley 2008, 2010, 2012). 
After summarizing several issues of growth at semimajor axes a ~ 125-250 AU (§2), we 
describe numerical simulations of planet formation for a variety of initial particle sizes and 
input parameters (§3-4). Within a ring of width 5a ~ 0.2a and total mass M 0 « 15 M®, 
ensembles of 1 cm to 10 m particles grow into super-Earth mass planets on time scales of 1- 
2 Gyr at 125-250 AU. After considering the implications of these results (§5), our discussion 
concludes with a brief summary (§6). 


2. BACKGROUND 

In the current core accretion theory, planets form in circumstellar disks of gas and 
dust (e.g., Goldreich et al. 2004; Chiang & Youdin 2010; Youdin & Kenyon 2013). After 
micron-sized grains within the disk agglomerate into cm-sized objects (e.g, Birnsticl et al. 
2010), continued coagulation (e.g., Windmark et al. 2012; Garaud et al. 2013) or some type 
of instability (e.g., Youdin & Goodman 2005; Johansen et al. 2012) concentrates solids into 
km-sized or larger planetesimals. When most of the solids end up in planetesimals, these 
objects collide, merge, and grow into planets (e.g., Kenyon & Bromley 2006; Chambers 2008; 
Kenyon & Bromley 2008; Kobayashi et al. 2010; Raymond et al. 2011). In a popular variant 
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of this picture, coagulation or instabilities produce only a few large planetesimals, which 
then rapidly accrete cm-sized pebbles (Ormel & Klahr 2010; Bromley & Kenyon 2011a; 
Lambrechts & Johansen 2012, 2014; Chambers 2014). 

Here, we focus on the growth of large swarms of planetesimals into super-Earth mass 
planets. When planetesimals form, their relative velocities are small. Physical collisions 
among planetesimals yield larger, merged objects. As planetesimals merge into larger and 
larger protoplanets, dynamical friction and viscous stirring excite the orbits of all the solids 
(e.g., Wetherill & Stewart 1989; Ida & Makino 1992a,b; Wetherill & Stewart 1993; Goldreich 
et al. 2004). Relative collision velocities grow in step with the escape velocity of the largest 
protoplanet. Eventually, collisions between small planetesimals produce debris instead of 
mergers. Collisions between debris particles make even more debris. The resulting collisional 
cascade grinds the debris into small dust grains; radiation pressure ejects these grains from 
the nascent planetary system (e.g., Dohnanyi 1969; Burns et al. 1979; Williams & Wetherill 
1994). Together, the collisional cascade and radiation pressure typically remove 80% to 95% 
of the initial mass in solids from the outer disk of a planetary system (Kenyon & Bromley 
2002 b, 2008, 2010; Kobayashi et al. 2010). As a result, the largest planets have maximum 
radii of ~ 3000-7000 km and maximum masses of 0.1-0.5 M e (e.g., Kenyon 2002; Kenyon 
& Bromley 2002b; Inaba et al. 2003; Kenyon & Bromley 2008, 2010; Kobayashi et al. 2010). 

Within this evolutionary sequence, there are several ways to limit the loss of small 
objects and to promote the formation of super-Earth mass planets at 125-250 AU. 

• Increase the binding energy Q* D of small planetesimals. Larger binding energies delay 
the collisional cascade, allowing protoplanets to grow to larger sizes. 

• Change the size distribution of the fragments produced in a collision. If the fragments 
have a steeper size distribution, it takes longer for the collisional cascade to grind them 
to dust. Longer collisional cascades enable larger protoplanets. 

• Include collisional damping of the smallest objects. If collisional damping can reduce 
the velocities of small collision fragments, collisional grinding proceeds more slowly 
and protoplanets grow larger. 

Of these options, modifying Q* D seems the least realistic. Modest changes to Q* D con¬ 
sistent with numerical simulations of the binding energy for small objects (e.g., Benz & 
Asphaug 1999; Leinhardt & Stewart 2009, 2012) lead to factor of ~ 2 changes in the mass of 
the largest object (e.g., Kenyon & Bromley 2010; Kobayashi et al. 2010; Kenyon & Bromley 
2012 ). Achieving much larger protoplanets requires factor of > 10 increases in the binding 
energy of 0.1-10 km objects. Calculations of the likely structure of these objects appear to 
preclude such large increases in Q* D . Thus, we eliminate this possibility. 
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Modifying the size distribution of the fragments is more promising. In all numerical 
simulations of planet formation, the size distribution of the debris is a power law n(r) oc r~ q 
with a maximum size r maxx i and a minimum size r minc ] (e.g., Wetherill & Stewart 1993; 
Kenyon & Luu 1999; O’Brien & Greenberg 2003; Kenyon & Bromley 2004a; Chambers 2008; 
Morbidelli et ah 2009; Kobayashi et ah 2010; Weidenschilling 2010). If m esc is the mass 
in fragments, defining r maXt d and q establishes the mass distribution uniquely (e.g., Kenyon 
et al. 2014). Many investigators set r max j ~ 7 r esc where r esc is the radius of a particle 
with mass m esc and 7 ~ 0.1-0.3. Because the equilibrium size distribution for a collisional 
cascade has q ~ 3.5-3.7 (Dohnanyi 1969; Williams & Wetherill 1994; O’Brien & Greenberg 
2003), many investigators fix q = 3.5. Although others employ more complicated formulae 
for r maXi d and q (Morbidelli et al. 2009; Weidenschilling 2010), these also require q < 4.0 at 
small sizes. 

In contrast to these restrictive approaches, detailed numerical simulations of binary 
collisions between planetesimals often yield debris with steep size distributions (e.g., Durda 
et al. 2004, 2007; Lcinhardt & Stewart 2012). For planetesimal radii Rq ^ 10-100 km, the 
debris has typical r max/ j ~ 0.1-1.0 R 0 . The smallest debris particles have sizes equal to the 
smallest size resolved in the calculation, r mtn/ j ~ 0.1-1 km. Over this range of sizes, the 
slope of the size distribution depends on the ratio r c of the collision energy to the binding 
energy. Collisions with larger r c have larger q. For r c > 1, extreme values for the slope of 
the size distribution are q ~ 4-6 (Leinhardt & Stewart 2012) or q m 8-10 (Durda et al. 2004; 
Morbidelli et al. 2009). 

In addition to exploring calculations with a range in q, it is important to consider the 
impact of collisional damping on the cascade (Goldreich et al. 2004). For typical surface 
densities of solid particles, collisional damping is effective for particle sizes of 1-10 mm 
(see also Kenyon & Bromley 2002a, 2014). Numerical simulations which follow such small 
particles generally do not include collisional damping (e.g., Stark & Kuchner 2009; Thebault 
et al. 2012 ; Gaspar et al. 2013). Calculations which include collisional damping rarely follow 
particles with sizes smaller than 1 cm (e.g., Wetherill & Stewart 1993; Kenyon & Bromley 
2008; Morbidelli et al. 2009; Kobayashi et al. 2010; Weidenschilling 2010; Kobayashi & Lohne 
2014). Thus, it is possible that a full evolutionary calculation with particle sizes smaller than 
1 cm will yield new phenomena. 

In the following section, we outline the tools we use to calculate the evolution of a swarm 
of particles with sizes of 1 jum to > 10 4 km. The Appendix summarizes several tests of our 
algorithms and comparisons with results of other approaches. Once the approach is set, §4 
describes results for evolutionary calculations at 125-250 AU. 
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3. PLANET FORMATION CALCULATIONS 

Over the past decade, we have developed Orchestra , an ensemble of computer codes for 
the formation and evolution of planetary systems. Orchestra consists of a radial diffusion 
code which derives the time evolution of a gaseous or a particulate disk (Bromley & Kenyon 
2011 a), a multiannulus coagulation code which infers the time evolution of a swarm of 
planetesimals (Kenyon & Bromley 2004a, 2008, 2012), and an n-body code which follows 
the orbits of gravitationally interacting protoplanets (Bromley & Kenyon 2006, 2011b, 2013). 
Within the coagulation code, Orchestra includes algorithms for treating interactions between 
small particles and the gaseous disk (e.g., Adachi et al. 1976; Weidenschilling 1977) and 
between coagulation mass bins and n-bodies (Bromley & Kenyon 2006). To treat interactions 
between small particles and the n-bodies more rigorously, the n-body calculations include 
tracer particles. Massless tracer particles allow us to calculate the evolution of the orbits 
of small particles in response to the motions of massive protoplanets. Massless and massive 
tracer particles enable calculations of the response of n-bodies to the changing gravitational 
potential of small particles (Bromley & Kenyon 2011a,b, 2013). 


3.1. Numerical Grid 

To explore a reasonable parameter space with good statistics, we perform coagulation 
calculations within a single annulus with semimajor axis a and width 5a around a star 
with mass M* = 1 M 0 . Within this annulus, there are M mass batches with characteristic 
mass rn k and radius r k (Wetherill & Stewart 1993; Kenyon & Luu 1998). Batches are 
logarithmically spaced in mass, with mass ratio 5 = m k+ i/m k . Each mass batch contains N k 
particles with total mass M k and average mass fh k = M k /N k . Particle numbers N k < 10 15 
are always integers. Throughout the calculation, the average mass is used to calculate 
the average physical radius f k , collision cross-section, collision energy, and other necessary 
physical variables. As mass is added and removed from each batch, the average mass changes 
(Wetherill & Stewart 1993). 

For any 5, numerical calculations lag the result of an ideal calculation with infinite mass 
resolution (see the Appendix). At a = 125-250 AU, calculations with 5 = 1.05-1.19 yield 
better solutions to the growth of 100+ km objects than calculations with 5 = 1.41-2.00. 
Although simulations with 5 = 1.05-1.10 allow somewhat better tracking of the late stages 
of planet growth, the improvement over results with 5 = 1.19 does not compensate for the 
factor of 2-4 extra cpu time required to complete these calculations. Thus, we consider a 
suite of calculations with 5 = 1.19 (= 2 1 / 4 ). 

In this suite of calculations, we follow particles with sizes ranging from a minimum size 
Knm = 1 /mi to the maximum size in the annulus r max . The algorithm for assigning material 
to the mass bins extends the maximum size as needed to accommodate the largest particles. 
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Typically, the number of large objects with r max > 1000 km in each calculation is small 
N k < 10; each protoplanet generally lies in its own mass bin. To conserve cpu time, we do 
not promote these objects into the n-body code. 


3.2. Initial Conditions 

All calculations begin with a mono-disperse swarm of planetesimals with initial size ro 
and mass density p p = 1.5 g cm -3 . These particles have initial surface density E 0 , total 
mass M 0 , and horizontal and vertical velocities h 0 and ty relative to a circular orbit. The 
horizontal velocity is related to the orbital eccentricity, e = 1.6 ( h/Vx ) 2 , where Vk is the 
circular orbital velocity. The orbital inclination is sin i = \/2v/Vk- 

For conditions at 125-250 AU, typical growth times arc longer than the 1-5 Myr lifetime 
of the gaseous disk (Haisch et al. 2001; Williams & Cieza 2011; Cloutier et al. 2014). Thus, 
we set the initial surface density of the gas to zero and ignore gas drag on small solids 
(Adachi et al. 1976; Weidenschilling 1977; Rafikov 2004). Compared to calculations which 
include gas drag (e.g., Wetherill & Stewart 1993; Chambers 2008; Kenyon & Bromley 2008; 
Kobayashi et al. 2010), our simulations retain a larger fraction of small particles which might 
experience significant radial drift during the likely lifetime of the gaseous disk. However, our 
goal is to isolate initial conditions and physical processes which might allow formation of 
super-Earths at 125-250 AU. We plan to consider the impact of gas drag in future papers. 


3.3. Evolution 

The mass and velocity distributions of the planetesimals evolve in time due to inelastic 
collisions, drag forces, and gravitational encounters. As summarized in Kenyon & Bromley 
(2004a, 2008), we solve a coupled set of coagulation equations which treats the outcomes 
of mutual collisions between all particles in all mass bins. We adopt the particle-in-a-box 
algorithm, where the physical collision rate is navf g , n is the number density of objects, cr 
is the geometric cross-section, v is the relative velocity, and f g is the gravitational focusing 
factor (Wetherill & Stewart 1993; Kenyon & Luu 1998). Depending on physical conditions 
in the disk, we derive f g in the dispersion or the shear regime (Kenyon & Luu 1998; Kenyon 
& Bromley 2012). For a specific mass bin, the solutions include terms for (i) loss of mass 
from mergers with other objects and (ii) gain of mass from collisional debris and mergers of 
smaller objects. 

As discussed in Wetherill & Stewart (1993), the most massive protoplanets on roughly 
circular orbits are ‘isolated’ from one another. Isolated protoplanets can accrete smaller 
objects but cannot collide with other isolated protoplanets. To identify isolated protoplanets, 
we find an ensemble of the largest objects whose combined gravitational range is smaller than 
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the width of the annulus. For an individual protoplanet, the gravitational range is 

Rg,k h (i H11 j.'h 2 ne/j (1) 

Here, K = 2\/3 and Rhm = [(mfe-t-m*,)/3M 0 ] 1//3 is the mutual Hill radius. The gravitational 
ranges of the isolated protoplanets satisfy the summation, 

kmax 

^ ^ ^ kRg,k — 3d i ( 2 ) 

k 

™min 

where is the number of protoplanets in mass bin k. 

Collision outcomes depend on the ratio Q c /Q*di where Q* D is the collision energy needed 
to eject half the mass of a pair of colliding planetesimals to infinity and Q c is the center of 
mass collision energy (see also Wetherill & Stewart 1993; Williams & Wetherill 1994; Tanaka 
et al. 1996; Stern & Colwell 1997; Kenyon & Luu 1999; O’Brien & Greenberg 2003; Kobayashi 
& Tanaka 2010). Following Wetherill & Stewart (1993), two colliding planetesimals with 
horizontal velocity hi, h 2 and vertical velocity V\, v 2 have relative horizontal and vertical 
velocities h c = {h\ + h 2 2 ) 1 ^ 2 and v c = (v 2 + v f) 1 / 2 (see also Kenyon & Luu 1998). The escape 
velocity of the colliding pair is v esc = ( 2Gm c /r c ) x / 2 , where m c = rn i + m 2 is the combined 
mass and r c = r\ + r 2 is the combined radius. The center of mass collision energy is then 

Qc = 0-5 fi(h 2 c + v 2 + v 2 esc )/m c (3) 

where // = m\m 2 /m c is the reduced mass. 

Consistent with N-body simulations of collision outcomes (e.g., Benz & Asphaug 1999; 
Lcinhardt et al. 2008; Lcinhardt & Stewart 2009), we set 

Q*d = Qb rfib + QgPp^ 9 (4) 

where Qyr^ b is the bulk component of the binding energy, Q g p g r^ 9 is the gravity component 
of the binding energy, and r is the radius of a planetesimal. 

In this study, planetesimals have Qb = 2 x 10 5 erg g' 1 cm 0 ’ 4 , j3b = —0.40, Q g = 0.22 erg 
g~ 2 cm 1 ’ 7 , and f3 g = 1.30. These parameters are broadly consistent with published analytic 
and numerical simulations (e.g., Davis et al. 1985; Holsapplc 1994; Love & Ahrens 1996; 
Housen & Holsapplc 1999). At small sizes, they agree with results from laboratory (e.g., 
Ryan et al. 1999; Arakawa et al. 2002; Giblin et al. 2004; Burchcll et al. 2005) and numerical 
(e.g., Leinhardt & Stewart 2009) experiments of impacts between icy o bjects. For the 
weakest objects with r ~ 10-100 m, our adopted Q* D is a factor of 3-10 smaller than in other 
studies (e.g., Bottke et al. 2010). Thus, our small planetesimals are relatively weak and easy 
to break. 

For two colliding planetesimals with masses m\ and m 2 , the mass of the merged plan¬ 
etesimal is 


m = m i + m 2 — m esc , 


( 5 ) 
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where the mass of debris ejected in a collision is 

m e s C = 0.5 (mi + m 2 ) • ( 6 ) 

The exponent 6 is a constant of order unity (e.g., Davis et al. 1985; Wetherill & Stewart 
1993; Kenyon & Luu 1999; Benz & Asphaug 1999; O’Brien & Greenberg 2003; Lcinhardt & 
Stewart 2012). In previous calculations, we adopted b = 9/8; here, we also consider 6 = 1 
(see also Kobayashi & Tanaka 2010; Kobayashi et al. 2010). 

To place the debris in the grid of mass bins, we set the mass of the largest collision 
fragment as rn maXy( i = 0.2 m esc and adopt a differential size distribution n(r) oc r~ q . After 
placing a single object with mass rn max> d in the grid, we place material in successively smaller 
mass bins until (i) the mass is exhausted or (ii) mass is placed in the smallest mass bin. Any 
material left over is removed from the grid. To explore a broad range of q (§2), we derive 
results for q = 3.5, 4.0, 4.5, 4.5, 5.0, and 5.5. 

To illustrate the general impact of 6 and q on the evolution, we consider two simple 
examples. When Q c /Q*d = 0-5 (early in the evolution), our approach yields ( m esc , rn max j) 
= (0.25, 0.05) m c when 6 = 1 and (0.23, 0.045) m c when 6 = 9/8. The larger exponent results 
in less debris which is placed in bins with smaller average mass. When Q C /Q* D = 1.5 (late 
in the evolution), ( m esc , rri maxx j) = (0.75, 0.15) m c when 6=1 and (0.79, 0.16) m c when 6 = 
9/8. The larger exponent produces more debris which is placed in bins with larger average 
mass. More debris makes collisional damping and dynamical friction more effective, which 
enhances the growth of the largest objects. Thus, the smaller (larger) exponent tends to 
enhance (retard) growth early in the evolution at the expense of slower (more rapid) growth 
later in the evolution. 

The broad range of q we consider has a similar impact on the timing and magnitude of 
runaway growth. When the size distribution is shallow (q = 3.5-4.0), debris is placed in a 
broad range of mass bins, limiting the impact of collisional damping and dynamical friction. 
When q is larger, debris is placed in a narrower range of mass bins, enhancing the impact of 
collisional damping and dynamical friction. Overall, we then expect calculations with larger 
q and larger 6 to have faster growth; calculations with smaller q and smaller 6 have slower 
growth. In §4, we discuss the impact of these choices in more detail. 

As we place the debris in specific mass bins, we also redistribute the kinetic energy per 
unit mass of each colliding pair of planetesimals. As in Kenyon & Luu (1998), we assume all 
collisions between mass batches conserve the horizontal and vertical components of kinetic 
energy. For an initial kinetic energy, m\{h\ + vf) + m 2 {h\ + u|), any merged planetesimal 
with mass m receives a fraction m/{m\ -\-m 2 ) of this kinetic energy; any fragment with mass 
nif receives a fraction m//(mi + m 2 ). Recalling the center of mass collision energy from 
eq. (3), this approach assumes that the escape velocity component of the collision energy is 
equal to the energy required to disperse the fragments to infinity. 



10 


To compute the evolution of the velocity distribution, we also include collisional damping 
from inelastic collisions and gravitational interactions. For inelastic and clastic collisions, we 
follow the statistical, Fokker-Planck approaches of Ohtsuki (1999) and Ohtsuki et ah (2002), 
which treat pairwise interactions (e.g., dynamical friction and viscous stirring) between all 
objects. For evaluating these interactions within a single annulus, we eliminate terms to 
calculate the probability that objects in one annulus interact with objects in other annuli 
(Kenyon & Bromley 2001, 2004b, 2008). We also compute long-range stirring from distant 
oligarchs (Weidenschilling 1989). The Appendix describes several tests of these algorithms. 


4. EVOLUTION OF THE LARGEST OBJECTS 

To evolve a sea of mono-disperse planetesimals in a single annulus, we set the size ro, 
the surface density E 0 , and the orbital elements eo and % o. For these simulations, r 0 = 
10 n cm where n is an integer between 0 and 8 inclusive. Beyond the snow line, the standard 
minimum mass solar nebula has E « 30 g cm -2 (a/1 AU) -3 / 2 (Kenyon & Bromley 2008), 
which implies E = 0.0215 g cm -2 at a = 125 AU. In an annulus with 5a = 0.2a, the total 
mass is Mo = 15.8 M®. For simplicity, we adopt the same mass in an annulus at 250 AU. 
In this annulus, E = 5.37 x 10~ 3 g cm -2 , roughly 70% of the surface density of a minimum 
mass solar nebula extrapolated to 250 AU. Previous calculations suggest the growth time for 
large objects scales inversely with the mass in solid objects (Kenyon & Luu 1999; Kenyon & 
Bromley 2002b, 2008, 2010). Thus, we consider only one surface density. 

In most published simulations, the initial orbital elements of planetesimals are set to 
match the escape velocity of the largest objects (e.g., Kobayashi et al. 2010; Weidenschilling 
2010). For particles with r 0 > 10-100 m at 125-250 AU, the time scale to reach equilibrium 
is longer than the 1-5 Myr lifetime of the gaseous disk (Appendix). Here, we set % = eo/2 
and adopt eo = 10 -4 for ro = 1 cm to 1 km, eo = 10~ 3 for ro = 10 km, and eo = 10~ 2 for 
r 0 = 100-1000 km. For the adopted E at 125-250 AU and these eo, swarms of particles are 
gravitationally stable and have eo intermediate between the likely e for particles within a 
turbulent gaseous disk and an equilibrium e set by a balance between stirring and damping 
(Appendix). After describing results for these initial conditions, we consider how different 
choices impact the results. 

In the next two sub-sections, we focus on the evolution of the largest objects at 125 AU 
(§4.1) and at 250 AU (§4.2). We then conclude this section with a brief discussion of the 
limitations of our calculations. 
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4.1. Results at 125 AU 

4-1.1. Super-Earth Formation 

Without gas drag, all calculations follow a similar pattern (e.g., Kenyon 2002; Kenyon 
& Bromley 2002b, 2004a, 2008, 2010). Starting from a mono-disperse set of planetesimals, 
collisions produce mergers and negligible debris. Growth of larger planetesimals allows dy¬ 
namical friction to circularize the orbits of the largest particles at the expense of raising 
the e and i of the smallest particles. Collisional damping gradually reduces the velocities of 
the smallest objects. Eventually, significant gravitational focusing factors enable runaway 
growth, where the largest planetesimals grow rapidly by factors of 10-100. As these pro¬ 
toplanets grow, viscous stirring raises the velocities of smaller planetesimals. Gravitational 
focusing factors diminish considerably. Runaway growth ends; oligarchic growth begins. 
Throughout oligarchic growth, collisions produce more and more debris. Once protoplanets 
reach radii of 1000 km or larger, the collisional cascade tries to remove nearly all of the small 
particles remaining in the annulus. As the cascade proceeds, protoplanets reach a maximum 
radius which remains fixed for the rest of the calculation. 

When 6 = 9/8 and q = 3.5, protoplanets reach radii of 1000 km or larger in roughly 1 Gyr 
(Fig. 1). In all calculations, growth starts slowly. Ensembles of 1 cm particles take ~ 10 5 yr 
to reach 10 cm sizes and another ~ 10 6 yr to reach 1 m sizes. After a very brief period 
of runaway growth, the largest objects have r max ~ 0.1-0.3 km and make the transition to 
oligarchic growth. Large protoplanets then reach 100 km sizes in 200-300 Myr. During a 
second period of runaway growth at 300-500 Myr, protoplanets grow past 1000 km. Some 
reach sizes of 5000-7000 km. We address the origin of this second phase of runaway growth 
in §4.1.2. 

Calculations starting with larger planetesimals follow the same steps. For r 0 = 1 cm to 
1 km, the time scale to reach the first runaway growth phase increases roughly linearly with 
the initial particle size. When tq = 10-1000 km, the largest planetesimals endure a long 
phase of slow growth, make a smooth transition to oligarchic growth, and then gradually 
accumulate enough material to reach sizes of 500-3000 km. 

This evolution is fairly independent of b (Fig. 2). At the smallest sizes (r 0 ~ 1 cm to 
1 km), the evolution of r max (t ) for b = 1 is nearly identical to the evolution for b = 9/8. 
Because planetesimals with r^ = 100-1000 km hardly grow in 10 Gyr, this evolution is also 
insensitive of b. For 10 km objects, however, growth is somewhat faster when 6=1 than when 
6 = 9/8. In our suite of calculations, it takes 150-175 Myr (225-275 Myr) to produce 100 km 
objects and 1000-1250 Gyr (1250-1250 Myr) to produce 1000 km objects when 6 = 1 (6 = 
9/8). The relative abundance of lower mass fragments causes this difference. During slow 
growth, calculations with 6=1 produce somewhat more fragments than calculations with 6 
= 9/8. The larger population of fragments circularizes the orbits of the largest protoplanets 
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more effectively, enabling a stronger runaway growth phase. Despite this difference, the 
largest protoplanets reach sizes of roughly 5000 km for all 6. 

To test the sensitivity of these results to stochastic variations within each simulation, 
we perform 5-10 calculations for q = 3.5 and each combination of b and r 0 . For b = 9/8 and 
ro = 1 cm to 0.1 km, the largest objects have radii r max « 3000-4500 km, with a median size 
r max ~ 3500 km. With r 0 = 1 km, the median r max ~ 6900 km; the median r rnax drops to 
4700 km for r 0 = 10 km. Ensembles of larger planetesimals grow very little: r max = 500 km 
for ro = 100 km and r max = 1500 km for r 0 = 1000 km. 

When 6=1, protoplanets typically reach smaller maximum sizes. Among the suite 
of calculations, variations in r max are similar to those for b = 9/8, 8% to 10% for sets of 
5-10 calculations with identical starting conditions. However, the median r max is 10% to 
15% smaller. When 6 = 1, two factors conspire to produce smaller protoplanets. During 
runaway and the early stages of oligarchic growth, collisions produce more debris when b 
— 1. At this time, collisional damping is ineffective. Thus, calculations with 6=1 lose 
somewhat more mass from the grid than those with 6 = 9/8. Throughout the late stages 
of oligarchic growth and the second phase of runaway growth, relative particle velocities are 
much larger. Collisions then produce more debris when 6 = 9/8. With so much debris, 
collisional damping becomes more effective for calculations with 6 = 9/8 than those with 6 = 
1. Damping promotes the growth of super-Earths, which are therefore more common when 
6 = 9/8. 

Although protoplanets with r max = 3000-7000 km are characteristic of growth with q 
= 3.5, the largest protoplanets in several calculations reach sizes exceeding 10 4 km (m > 
1 M e ). Each of these simulations yields 2-3 super-Earths. For ro = 1 cm to 10 km, 3% (2%) 
of simulations with 6 = 9/8 (1) produce super-Earth mass planets. Each of these examples 
began with small planetesimals, ro < 100 cm; however, small number statistics precludes 
correlating occasional super-Earth formation with ro- 

Super-Earth formation is much more common in calculations with q > 4.5 (Fig. 3). For 
all q, protoplanets grow at similar rates for 100 Myr to 1 Gyr. At later times, protoplanets 
in calculations with q > 4.5 grow much more rapidly. In simulations with r 0 = 1 cm to 
1 km and 6 = 9/8, protoplanets reach median maximum sizes r max « 9000-12000 km (m ~ 
0.75-1.8 The variations about the medians are much larger, ~ 20% to 25% instead of 
the 10% range for q = 3.5-4.0. When r 0 = 10 km, maximum sizes are roughly 20% smaller. 
None of the simulations with larger planetesimals (r 0 > 10 km) yields a super-Earth. 

Results for 6=1 are similar. Suites of simulations with identical starting conditions 
have typical protoplanet radii of 7500-10000 km, roughly 15% smaller than calculations with 
6 = 9/8. Because the variations about the median size are large, ~ 25%, many simulations 
with 6 = 1 and r 0 = 1 cm to 1 km produce super-Earth mass planets. Ensembles of larger 
protoplanets with r 0 > 10 km never produce super-Earths. 
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Figs. 4-5 summarize results for r max as a function of 5, ro, and q. When ro = 1 cm to 
0.1 km, b = 9/8, and q = 3.5 or 4, fewer than 5% of all simulations produce super-Earth 
mass planets (Fig. 4). Although rings composed of larger planetesimals (r 0 = 1 km) usually 
yield more massive planets when q = 3.5 or 4, super-Earths are still rare. In contrast, super- 
Earths form in more than half of the simulations with ro = 1 cm to 1 km and q = 4.5, 5, 
or 5.5. When 6 = 1 and q = 4.5, 5, or 5.5, super-Earth formation is somewhat less common 
(Fig. 5). With r 0 < 1 km, super-Earths form in roughly 25% of the simulations with q > 
4.5. Among simulations with q < 4, super-Earth formation is rare. 


4-1-2. Collisional Damping when q = 3.5-4-0 


Compared to previous calculations of planet formation at 100-150 AU with similar 
starting conditions and fragmentation parameters, our new calculations yield fairly similar 
results. Published multiannulus calculations at 100-150 AU with identical fragmentation 
parameters and q = 3.4 have a median r max fa 1250 km (Kenyon & Bromley 2010). In 
the calculations for this paper with q = 3.4 (3.5), the median r max is 1500-2000 km (2000- 
3000 km). By neglecting gas drag and including the evolution of small particles with radii 
smaller than 1 m, the calculations in this paper produce somewhat more massive large 
objects. 


Although the new calculations allow the growth of larger objects when q = 3.4-3.5, the 
collisional cascade still converts most of the initial solid mass into particles with sizes smaller 
than 1 /jrn. Radiation pressure rapidly ejects these particles. In Kenyon & Bromley (2010), 
the cascade removes from 0.1% (r 0 = 100 km) to 90% (r 0 = 1 km) of the initial mass. Our 
new results suggest a weak dependence of mass loss on the initial planetesimal radius for 
small planetesimals; the cascade removes 80% of the initial mass when ro = 1 cm and more 
than 90% when r o = 1 km. When ro > 10 km, the cascade removes very little mass as in the 
Kenyon & Bromley (2010) simulations; the slow evolution of swarms of large planetesimals 
then prevents the formation of large planets (see also Kenyon & Bromley 2010). 


Despite the rough similarity in r max and the ejected mass among these simulations, the 
time evolution of the size distribution is very different. In previous multiannulus simulations, 
the cumulative size distribution can be split into three power-laws (e.g., Kenyon & Bromley 
2004d): 


r 

n c (r) = < 




n s r~ qs 

r < r 2 

ni'r~ qi 

ry < r < ri 

ni'r~ qi 

r > r i 


(7) 


Here, r\ ~ 10-100 km, r 2 ~ 0.1-1 km, q s fa 2.5-2.7, g* ~ 0-2, and qi ~ 3-4 (see also Kenyon 
& Bromley 2012). Growth sets the slope at large sizes, r > rp collisional debris is responsible 
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for the slope at small sizes, r < r 2 (Kenyon & Bromley 2004d, 2008, 2012; Pan & Sari 2005; 
Schlichting & Sari 2011; Schlichting et al. 2013). At intermediate sizes, competition between 
collisional destruction and debris production sets the slope. 


To visualize the evolution of the size distribution in the new simulations, we normalize 
the cumulative size distribution to the expected size distribution at small sizes, n c;re i{r) = 
n c (r)fn s r~ qa . With this normalization, the expected cumulative size distribution is: 


1 

r < r 2 

n 2 r qs ~ qi 

r 2 < r < r\ 

ni r q s- q i 

r > r\ 


For this expression, we set ri\ = ni/n s and n 2 = rii/n s . 


( 8 ) 


Fig. 6 illustrates some of the advantages of the relative size distribution. In the lower 
panel, three cumulative size distributions with different power-law components look fairly 
similar. Shifting to the relative size distribution in the upper panel, the three curves are 
very distinct. With q t < q s < qi , the relative cumulative size distribution rises from r 2 to 
r\ and then falls from r\ to r max . This normalization makes it fairly easy to infer T\ and r 2 
from the breaks in the relative size distribution. 


Fig. 7 shows the evolution of the cumulative relative size distribution for the first 300 Myr 
of a calculation with ro = 1 cm, b = 9/8, and q = 3.5. At 0.1-1.0 Myr, the size distribution 
approximately follows the expected shape, with (i) a flat debris tail at r « 1-100 /mi; (ii) 
a sharp rise at r « 100 /mi to r ~ 1-10 cm; and (iii) a steep fall for r > 1-10 cm. The 
small dip at r ~ 100 /mi establishes particle sizes where removal of particles by merger with 
much larger particles dominates debris production by destructive collisions. The peak of the 
curves corresponds to particle sizes which contain most of the mass. 

As the calculation proceeds, larger and larger particles contain more and more of the 
mass. It takes ~ 30 Myr for the peak in the relative size distribution to reach 10 m and 
another 270 Myr to reach 100 m. The growth of the largest particles follows this evolution, 
with 10 m particles in 1 Myr, 1 km particles in 3 Myr, and 10 km particles in 30 Myr. By 
300 Myr, the largest particles have reached sizes of 200 km. 

Throughout this evolution, the size distribution develops several distinctive features. At 
small sizes (r « 10-100 /mi to 1-10 cm), the debris tail contains more particles and evolves 
towards a single power-law. Among larger particles with r « 1-100 km, there is a steeper 
power-law composed of large objects which grow by accreting much smaller objects (see also 
Kenyon 2002; Goldreich et al. 2004; Kenyon & Bromley 2008, 2010, 2012; Schlichting & Sari 
2011 ). 

Once the largest particles have sizes exceeding 500-1000 km, the size distribution 
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changes dramatically (Fig. 8). From 300 Myr to 500 Myr, (i) the clip at sizes of 10 m 
shifts to smaller and smaller particles and becomes progressively deeper; (ii) the peak at 
sizes of 100 m shifts to larger and larger sizes; and (iii) the fraction of mass in the debris 
tail grows significantly. These changes herald a short (second) runaway growth phase, where 
the largest particles with r « 100 km to 1000 km begin to contain most of the mass and 
achieve a size distribution with a standard slope qi ~ 3 (see also Kenyon & Bromley 2012). 
By rsj 1.5 Gyr, collisions have removed most of the mass in small particles, effectively ending 
growth of the largest particles. 

To understand the interesting evolution of the size distribution, it is helpful to con¬ 
sider the behavior of the velocity distribution. To simplify this discussion, we calculate the 
eccentricity e*, of each mass bin relative to the Hill eccentricity = (m max /3M Q Y^ 3 , 

(•k.rfl ('A -/('ll j (9) 

where m max is the mass of the largest object in the annulus. 

During the first period of runaway growth and oligarchic growth, the relative eccentricity 
distribution follows a standard pattern (Fig. 9; see also Wetherill & Stewart 1993; Kenyon 
& Luu 1999; Ohtsuki et al. 2002; Kenyon & Bromley 2002b; Goldreich et al. 2004; Kenyon 
& Bromley 2004a, 2008). When the largest objects reach sizes of 1 km, collisional damping 
and viscous stirring have established an equilibrium where the random velocities of small 
particles (r <0.1 km) are ~ ten times the Hill velocity of the largest particles, a factor of 2-4 
smaller than in our previous calculations where r m j n = 1 nr (e.g., Kenyon & Bromley 2008, 
2010 , 2012). At the largest sizes, dynamical friction reduces the eccentricities to a minimum 
£k,rei ~ 0.1-0.2. In between, dynamical friction maintains an approximate power-law, with 
?k,rei r~ n and n ~ 0.5-1. Despite the larger stirring as r max reaches 30 km, a rough 
balance between collisional damping and viscous stirring preserves (i) the constant relative 
eccentricity, e^rei = ek/eH ~ 10, among the small particles and (ii) the power law decline 
among the larger particles. 

As the evolution proceeds past 100 Myr, it is convenient to consider the behavior of 
large, intermediate, and small particles. For our adopted Q* D , large objects with r > 10 km 
do not shatter; they gain mass and energy from collisions with smaller particles in the grid. 
Growing protoplanets also try to stir all other particles to larger and larger velocities. The 
magnitude of the stirring grows with the size of the protoplanet; the impact on smaller 
particles is fairly independent of the size of the small particles. Dynamical friction with 
small objects circularizes the orbits of the protoplanets. 

Intermediate mass particles with r ~ 10-100 m contain most of the mass but are easily 
shattered. Collisions among these objects place mass and energy into smaller mass bins. 
Compared to stirring by large particles, damping is relatively inefficient. Thus, intermediate 
mass bins lose mass and gain velocity with time. 
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Small objects with r < 1 cm contain little mass but do not shatter. These particles gain 
mass from the debris of collisions among intermediate mass particles and lose mass from 
collisions with large objects. Slow growth among these particles gradually shifts material to 
larger bins. As this material moves into the intermediate mass bins, destructive collisions 
cycle it back into smaller mass bins. With most of the mass in intermediate mass particles, 
the small particles gain mass with time. 

The evolution of ek^ei for small particles depends on the relative importance of collisional 
damping and stirring by large particles. Without any gain in mass from destructive collisions, 
stirring by large particles dominates damping. However, destructive collisions rapidly cycle 
material from 10-100 m particles into smaller particles. Collisions among stronger small 
particles are not destructive. Thus, the net mass in small particles grows (Figs. 7-8). During 
this evolution, collisional damping gradually overcomes stirring by larger particles; 6k, re i 
gradually declines. 

Ultimately, the shape of the Q* D relation drives this evolution. In alternative calculations 
where Q* D oc r~^ b and [3b < 0 at small radii, dramatic reductions in ek, r ei among small particles 
cannot occur. When the largest particles reach 10-100 km sizes, the smallest particles in the 
grid shatter first. As large objects continue to grow, shattering impacts larger and larger 
particles, but the mass in the smallest particles always decreases. This evolution severely 
limits collisional damping and allows a collisional cascade to destroy all objects below some 
size limit. 

With (3b > 0, intermediate mass particles always shatter first. For our adopted pa¬ 
rameters, particles with r ~ 30 m have the smallest Q* D . As the random velocities of all 
particles slowly increase, debris from collisions of these weak planetesimals populates mass 
bins where particles are stronger and do not shatter. Collisional damping can then reduce 
£k,rei f° r objects in these mass bins. 

The drop in 6k,rei from collisional damping is responsible for a second runaway in the 
growth of the largest objects. With ek,rei ~ 0.5-2 at 300-500 Myr, gravitational focusing 
factors rise dramatically. As the largest protoplanets grow from 500 km to 2000 km, viscous 
stirring overcomes collisional damping for smaller and smaller objects. Thus, the ‘collisional 
damping front’ moves from roughly 10 m at 300 Myr to 10 cm at 500 Myr to 0.1 cm at 
1.5 Gyr; the relative eccentricities for the larger objects rise from ek, r ei ~ 5-6 at 500 Myr to 
€k,rei ~ 20 at 1.5 Gyr to ek^ei ~ 30-40 at 10 Gyr. 

In all simulations with q = 3.5-4.0, the ultimate sizes of the largest objects depend on 
the timing of the dramatic reduction in ek,rei for the smallest objects. Usually, oligarchic 
growth produces large objects with ek, r ei > 3-10 before collisional damping is able to reduce 
the relative eccentricities for the smallest objects. Collisions between small objects and 
oligarchs then occur in the dispersion regime, where gravitational focusing factors are small 
and growth is slow. Despite damping, collisions among small objects are destructive. The 
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cascade then slowly removes the mass in small objects before oligarchs reach super-Earth 
sizes. 

Sometimes, conditions conspire to enable the production of several super-Earth mass 
planets. When oligarchs and small objects have e^^ei 1-3, collisions are in the shear regime 
where gravitational focusing factors are large. Growth is then very rapid. For simulations 
with q = 3.5-4.0, the debris from every destructive collision is spread out over several orders of 
magnitude in radius. The large radius spread in the debris (i) limits the impact of dynamical 
friction on the oligarchs (enabling higher orbital eccentricities) and (ii) reduces the ability 
of collisional damping. Thus, super-Earth formation is very rare. 


4-1.3. Collisional Damping when q = 4-5-5.5 

When q > 4.5, the second phase of runaway growth becomes more explosive (Figs. 3). 
At 300 Myr to 1 Gyr, the largest particles rapidly grow from r max = 2000-4000 km to super- 
Earth sizes with r max > 10 4 km. Although super-Earth production is fairly insensitive to q 
for q > 4.5, maximum planet sizes depend on 6. Calculations with b = 9/8 almost always 
yield larger planets than calculations with 6=1. 

Two factors enable routine super-Earth formation when q > 4.5. In all calculations, 
the collisional cascade begins after the formation of large objects with r max > 500-1000 km. 
Debris production grows dramatically. Because radiation pressure ejects particles smaller 
than 1 /jrn, the mass loss rate from the grid also grows. Although the debris production rate 
is insensitive to q, calculations with larger q concentrate the debris into a smaller range of 
radii and eject fewer particles with radii smaller than 1 fin i. At similar stages of evolution, 
more mass remains in the grid. Thus, the largest protoplanets tend to grow to larger masses. 

Calculations with larger q also enhance the impact of collisional damping. When q 
is small (3.5-4.0), debris is spread out over a broad range of sizes. Collisional damping - 
which depends on the collision rate - then reduces the eccentricities of small particles more 
efficiently than those of larger particles (Fig. 10). With less mass in smaller objects than 
in larger objects, the growth of the largest protoplanets is fairly slow. When q is large 
(4.5-5.5), debris is concentrated in the largest particles. Collisional damping is then more 
effective for larger particles, which now contain much more mass than the smaller particles. 
Once relative eccentricities fall below the Hill velocity, protoplanets rapidly accrete the more 
massive particles and reach super-Earth sizes. 

Fig. 11 compares the relative size distributions for calculations with various q. To il¬ 
lustrate the differences clearly, we focus on evolution times where the largest objects have 
r max = 2000-3000 km and are just starting their second phase of runaway growth. When q 
= 3.5, the slope of the size distribution for small particles with r < 1-10 cm is shallow. Cal¬ 
culations with larger q have steeper size distributions with more mass in the largest objects. 
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Among larger particles with r > 10 4 cm, the size distributions are nearly independent of q. 
At intermediate sizes (10 cm < r < 10-100 m), the concentration of mass in larger particles 
is striking. 

Fig. 12 compares the relative eccentricity distributions for the same epochs as Fig. 11. 
At small sizes (r < 1-10 cm), the changing impact of collisional damping with increasing q 
is apparent: damping is more effective at small sizes when q = 3.5 and at large sizes when 
q = 4.5 and 5.5. With more debris in larger particles, the damping front - which marks the 
transition between particles with large e and small e - extends to larger r for calculations 
with larger q. For all particles with r > 1-10 cm, calculations with larger q have smaller 
relative eccentricities. When q = 3.5, the largest protoplanets have ek, re i > 1; despite the 
small ek t rei f° r the smallest particles, these planets grow relatively slowly during the second 
runaway growth phase. When q = 4.5-4.5, the largest protoplanets have e^ re i < 1; these 
planets can grow rapidly. 


4-1-4- Evolution with Different Starting Conditions 

To examine how the results of these simulations depend on initial conditions, we consider 
suites of simulations where the initial eccentricity and inclination are a factor of three larger 
or smaller than our nominal eo and i 0 . For each r 0 , we perform another 3-5 calculations 
with the different starting values of e and i. Within each suite of simulations, the median 
r max is nearly identical; there is no obvious variation of r max with eo- The range in r max for 
these simulations is also similar, ~ 20% to 25%. 

Within the full suite of simulations, there are modest differences in the time evolution 
of r max and the size distribution as a function of starting e and i. In rings composed of small 
objects with ro < 10-100 m, gravitational focusing is negligible at the start of the simulation. 
Thus, growth rates depend only on the geometric cross-section and the relative numbers of 
small and large objects. Although collisions produce more debris when the initial e and i 
are large, the debris contains a negligible amount of mass. Thus, the initial slow growth of 
small particles is independent of eo- 

As the evolution proceeds, the rapid growth of the largest particles has two regimes. 
When the largest objects make the transition from slow growth to runaway growth, gravi¬ 
tational focusing factors increase dramatically. Ensembles of particles with smaller initial e 
reach this transition earlier than ensembles of particles with larger initial e. Defining tn as 
the time of this transition, our results suggest tn(e o/3) ~ 0.25f#(3eo). 

When the largest objects make the transition from runaway to oligarchic growth, the 
growth rate again depends on the initial eccentricity. Throughout runaway growth, collisions 
in simulations with the largest eo produce the most debris. More debris leads to more effective 
collisional damping and somewhat more effective dynamical friction between the largest and 
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smallest objects. As a result, oligarchs grow faster in calculations with larger eo- Despite 
starting out more slowly, these oligarchs overtake the oligarchs in calculations with smaller 
eo- Thus, calculations with larger eo produce super-Earths somewhat more rapidly than 
calculations with smaller eo- 

Among swarms of large particles with ro > 1 km, the evolution seems to depend little 
on eo• Over typical growth times of 10-100 Myr for 1 km particles, collisional damping, 
dynamical friction, and viscous stirring produce similar size distributions for calculations 
with different eo- At the start of runaway growth, all calculations then appear to follow 
a similar sequence. Within this sequence, there is a large dispersion of outcomes due to 
stochastic variations in the collision rate. This dispersion dominates small differences in 
outcomes due to the starting conditions. 

Fig. 13 illustrates some of these points for calculations with 6 = 1 and q = 4.5. When 
ro = 10 cm, particles with smaller initial e grow somewhat faster. Once r max exceeds 10 km, 
oligarchs in calculations with larger e grow faster and eventually overtake oligarchs with 
smaller initial e. At late times, the largest objects enter a second runaway growth phase, 
where they accrete material rapidly at rates independent of e 0 . Eventually, the largest objects 
reach similar (super-Earth) sizes. When ro = 1 km, growth is more stochastic. Sometimes, 
small e calculations lag those with larger e (as in the Figure). In other simulations, large e 
calculations lag. All eventually produce 1-2 protoplanets with large radii. 


4.2. Results at 250 AU 

At 250 AU, the evolution from small planetesimals into snper-Earth mass planets follows 
the behavior described for calculations at 125 AU. To avoid repeating the discussion of §4.1, 
we summarize differences between the two sets of simulations in the rest of this section. 

The time scale for planet formation is the simplest difference between sets of calculations 
at 125 AU and at 250 AU. At any semimajor axis, the time scale for planets to reach 
identical sizes scales inversely with the initial mass in planetesimals and linearly with the 
orbital period. With identical initial planetesimal masses in our two sets of calculations, the 
formation time scales with orbital period. Thus, planets at 250 AU take roughly 3 times 
longer to form than those at 125 AU. For calculations with r 0 = 1 cm, this difference is 
not a major issue: 1 Gyr formation times at 125 AU become 3 Gyr time scales at 250 AU. 
Both time scales are shorter than the age of the Solar System. For calculations with ro 
~ 10-100 m, this difference is crucial: 2-3 Gyr formation time scales at 125 AU become 
6-10 Gyr time scales at 250 AU. Thus, large particles at 250 AU do not form super-Earth 
mass planets over the age of the Solar System. 

Fig. 14 illustrates this point for the growth of 1 cm and 10 m particles. When ro = 
1-10 cm, growth follows the same stages outlined in §4.1. After 1-10 Myr of slow growth, 
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periods of runaway growth, oligarchic growth, and runaway growth produce super-Earth 
mass planets in 2-4 Gyr. The growth time is roughly 3 times longer than the time scale at 
125 AU (see Figs. 1-2). When r 0 > 10 m, every phase of growth takes longer. Thus, the 
first runaway growth phase begins at 100 Myr to 1 Gyr, instead of 1-10 Myr. Although 
collisional damping still promotes rapid growth during the second runaway growth phase, 
rapid growth begins at 5-10 Gyr - too late to allow super-Earth formation within the age 
of the Solar System. 

Although growth at 250 AU takes longer, small particles reach super-Earth sizes inde¬ 
pendently of q (Fig. 14). For ro = 1-10 cm and all g, collisional damping drives the second 
period of runaway growth, where 2000-3000 km protoplanets rapidly grow into super-Earths. 
The time scale to produce super-Earths is also insensitive to q, 2-4 Gyr. Although super- 
Earths form for all q, calculations with q = 4.5-5.5 typically yield planets roughly twice as 
massive as calculations with q = 3.5-4.0. 

At larger sizes (ro > 1-10 m), final planet masses depend on q. Calculations with q = 
3.5-4.5 stall at typical sizes of 4000-5000 km; simulations with q = 4.5-5.5 reach radii of 
7000-8000 km at 10 Gyr. On longer time scales, calculations with large q eventually reach 
super-Earth masses. Models with smaller q never produce massive planets. 

Results at 250 AU seem fairly insensitive to b. For ro = 1-10 cm, all calculations yield 
super-Earths. The largest planets have masses of 2-3 M e for 6=1 and b = 9/8. For 
r 0 > 10 m, planets also have maximum sizes that do not depend on b. When r 0 = 1 m, 
calculations with b = 9/8 and q = 4.5-5.5 routinely yield super-Earth mass planets on time 
scales of 2-4 Gyr. Simulations with b = 9/8 (1) and q = 3.5-4.0 (3.5-5.5) cannot produce 
super-Earths in 10 Gyr. 

Figs. 15-16 summarize r max for the complete suite of calculations at 250 AU. When b 
= 1, nearly all calculations with ro = 1-10 cm yields a super-Earth mass planet. Although 
roughly half of simulations with ro = 1 m yield a super-Earth, many produce smaller planets 
with ro = 3000-7000 km. These results are fairly insensitive to q. With ro = 10 m, only a 
few calculations yield super-Earths. Ensembles of larger planetesimals never produce super- 
Earths. When b = 9/8 and ro = 1 cm to 1 m, super-Earth formation is also common. Larger 
planetesimals almost never produce a super-Earth. These results are also fairly independent 
of q. 


4.3. Limitations of the Models 

In previous papers, we have summarized various limitations in coagulation calculations 
for planet formation (e.g., Kenyon & Luu 1999; Kenyon & Bromley 2004a, 2008). Here, we 
focus on how choices in methods and initial conditions impact our results. 
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4-3.1. Resolution 

Finite resolution - in the mass spacing parameter 5 and the width of the annulus 5a 
- is characteristic of all numerical coagulation calculations. As noted in the Appendix, 
simulations with finite 5 lag the results of a calculation with ‘perfect’ mass resolution (see 
also Wetherill 1990; Kenyon & Luu 1998). For 5 = 1.2, this lag is negligible relative to 
simulations with 5 = 1.05-1.10 (Appendix). Thus, uncertainties due to finite 5 are small. 

Adopting a single annulus for our calculations limits our ability to follow the evolution of 
a planetesimal swarm as a function of semimajor axis (e.g., Spaute et al. 1991; Weidenschilling 
et al. 1997; Kenyon & Bromley 2004a, 2008). Despite this disadvantage, single annulus 
calculations yield fairly accurate estimates for the growth times and outcomes of coagulation 
calculations. For the initial conditions adopted here, the formation time depends only on 
the orbital period. Across the width of our single annulus, 5a/a = 0.2, the orbital period 
changes by 35%. If the surface density is constant across the annulus, planetesimals at the 
inner edge of the annulus grow 35% faster than planetesimals at the outer edge of the grid. 
With typical formation times of ~ 1 Gyr and a factor of two range in the formation time, 
uncertainties due to the size of the annulus are small. Thus, our results yield good first 
estimates for the growth of super-Earths at 125-250 AU. 

Although we allow any maximum size in the coagulation grid, setting a hnite lower limit 
on particle size r mm changes the evolution of larger particles. During slow and runaway 
growth, particles with r ~ r mm have negligible mass and little impact on the evolution of 
larger particles. Once the collisional cascade begins, debris rapidly fills the smallest mass 
bins. Particles with sizes r ~ 10-100 r rnm are swept up by much larger protoplanets and 
destroyed by particles with somewhat smaller sizes r ~ 1-10 r min . With no particles smaller 
than r min in the grid, our calculations have fewer destructive collisions among particles with 
r % T m i n compared to a calculation with an infinitesimally small r m i n . Thus, our calculations 
produce an ‘excess’ of small particles relative to an ‘ideal’ calculation (see also Carnpo 
Bagatin et al. 1994; O’Brien & Greenberg 2003; Krivov et al. 2006, 2008, and references 
therein). 

The standard method to address this issue is to extend the size distribution below r m in, 
assume these particles have relative velocities similar to particles with r r mm , derive 
collision rates and outcomes with particles in the grid, and correct the excess number of 
particles with r cs r min (e.g., Carnpo Bagatin et al. 1994; O’Brien & Greenberg 2003). 
Although this approach is attractive when r m ; n ~ 1-100 cm, it is problematic when r mm 
fti 1 /jrn. For solar-type stars, radiation pressure likely places smaller particles on highly 
eccentric or hyperbolic orbits as long as the ensemble of solids is optically thin (e.g., Burns 
et al. 1979). Relative collision velocities are then much larger, but collision rates are much 
smaller. 

We plan to perform a set of multiannulus calculations to quantify outcomes with realistic 
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orbital geometry for small particles. To provide an initial estimate of the importance of the 
small size cut-off, we perform several calculations with r min = 0.1 /jrn. Aside from shifting 
the excess of small particles to smaller sizes, this change has a modest impact on the overall 
size distribution or the growth of small objects. Thus, this uncertainty is minimal. 


4-3.2. Fragmentatiori 

The fragmentation algorithm is another uncertainty. We use a standard energy scaling 
approach which mimics the results of more detailed SPH and n-body calculations of binary 
collisions. Although our fragmentation parameters - Qb, Q g , (3b, and (3 g - are consistent 
with the detailed calculations, other choices are possible (e.g., Morbidelli et al. 2009; Bottke 
et al. 2010; Kenyon & Bromley 2010; Kobayashi & Lohne 2014). Based on published results, 
it is straightforward to assess how alternative choices impact our results. For simplicity, 
we separate the discussion into the gravity regime for large objects (r >0.1 km) and the 
strength regime for small objects (r <0.1 km). 

In the gravity regime, there is general agreement among n-body and SPH calculations for 
/ 3 g (Benz & Asphaug 1999; Leinhardt & Stewart 2012). Small variations about the adopted 
value have little impact on collision outcomes. Current simulations suggest at least a factor of 
2-3 uncertainty in Q g (Morbidelli et al. 2009; Bottke et al. 2010). Although the ultimate sizes 
of protoplanets grow with increasing Q g (e.g., Kenyon & Bromley 2010), efficient collisional 
damping probably limits the importance of this parameter on our conclusions. Reducing Q g 
by a factor of 100 or more probably prevents the formation of 1000 km or larger protoplanets. 
More modest changes, however, probably have little impact on our results. 

In the strength regime, different choices for Qb and (3b do not affect runaway and oli¬ 
garchic growth (Kenyon & Bromley 2010). Making small particles harder to break (by in¬ 
creasing Q b and decreasing (3b) increases the population of small objects and the effectiveness 
of collisional damping, enabling a broader set of initial conditions to produce super-Earths. 

For smaller particles, likely outcomes depend on the relative importance of shattering 
and collisional damping. As an example, we consider the impact of reducing Q* D by a factor of 
ten for particles with r < 100 m in the strength regime. In our nominal calculations, collisions 
with impact velocities v > 400 cm s -1 shatter 30 m particles. As the largest particles grow 
past 100 km, collisional damping maintains collision velocities for these particles just above 
this limit (Fig. 9). Although these particles shatter, particles much smaller than 10 m and 
much larger than 100 m do not shatter, which allows collisional damping to initiate the 
second phase of runaway growth. When Q* D is a factor of ten smaller, impact velocities 
v > 125 cm s -1 shatter small particles, distributing debris into smaller mass bins. If the 
collisional damping rate remains identical to that in our nominal calculations, continued 
shattering removes the debris rapidly. Collisions of much more massive particles cannot 
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produce debris fast enough to replenish material lost from smaller mass bins. The collisional 
cascade then stalls the growth of the largest objects, which never reach super-Earth masses. 

ffowever, increased debris production associated with a smaller Q* D also increases the 
effectiveness of collisional damping. A simple estimate suggests that damping counters the 
impact of a larger rate of shattering. When Q* D is a factor of ten smaller, debris production 
for fixed Q c is a factor of ten larger (eq. 6). With a factor of ten more mass in small 
particles, the damping rate is a factor of ten larger. Setting the approximate rates for viscous 
stirring and dynamical friction from large particles to this damping rate (e.g., Goldreich et al. 
2004), the likely equilibrium velocity of small particles, 75-100 cm s -1 , is smaller than the 
125 cm s _1 shattering velocity. Thus, collisional damping appears capable of countering any 
large reduction in Q* D and maintaining the evolution on the path to super-Earth formation. 

As noted in §4.1.2, setting [5 g < 0 limits the impact of collisional damping on the velocity 
evolution of small particles. When /3 g < 0, the smallest particles in the grid shatter first. 
Debris from these collisions permanently leaves the grid. As successively larger particles 
begin to shatter, smaller particles continue to shatter. More debris leaves the grid. With too 
little mass among smaller particles, collisional damping is ineffective. Velocities for small 
particles remain large, preventing a second runaway growth phase leading to the formation of 
super-Earth mass planets. Maximum protoplanet radii are then small, ~ 2000 km. Repeating 
our calculations with different parameters for Q* D is necessary to confirm this conclusion. 

Aside from choosing q and the four parameters in the Q* D relation, we also set the 
mass of the largest fragment in the debris, rri^ = 0.2m esc , where m esc is the total mass of 
the debris (see also Wetherill & Stewart 1993; Kenyon & Luu 1999). This approach differs 
from Morbidelli et al. (2009), nij j a; 0,004e^ c / 4< %) 2 m esc ; Kobayashi & Tanaka (2010), ?tll ~ 
€ (Qc/Q*D) lm esc an d e ~ 0.01-0.5; and Weidenschilling (2010), rrii = 0.5 (Q c /Q* D )^°' 7S m esc . 
For moderate collision energies (Q c ~ Q* D ), our algorithm places a somewhat smaller fraction 
(Weidenschilling 2010), a similar fraction for e = 0.2 (Kobayashi & Tanaka 2010), or a much 
larger fraction (Morbidelli et al. 2009) of the ejected mass in the largest particle within the 
debris. At low collision energies ( Q c <C Q * D ), Kobayashi & Tanaka (2010) and Weidenschilling 
(2010) have more massive fragments; Morbidelli et al. (2009) has less massive fragments. At 
high collision energies ( Q c 3> Q * D ), all other algorithms place more material in smaller 
fragments. 

Although accurately assessing the impact of these choices requires a broad set of addi¬ 
tional simulations, we can outline the probable role of rriL in our calculations at 125-250 AU. 
The Morbidelli et al. (2009) algorithm always places more mass in less massive particles 
within the debris. Per collision, this approach thus assigns more mass to particles with 
r < r m i n , leading to less mass among debris with r > r m i n . Fewer debris particles results 
in less efficient collisional damping and limited growth during any second phase of runaway 
growth. Super-Earth formation seems unlikely. 
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In the Kobayashi & Tanaka (2010) and Weidenschilling (2010) approaches, super-Earth 
formation probably follows a similar path as in our calculations. At the start of the second 
runaway growth phase, most of the mass is in particles with r > 1 km and Q c ~ Q* D . For 
identical q, these algorithms place the debris in similar mass bins. Among small objects (r < 
1 cm) with Q c < Q* D , our algorithm removes mass from the grid more rapidly than these two 
approaches and thus underestimates the growth rate for super-Earths. At intermediate sizes 
(10 cm to 0.1 km), our approach places mass in larger fragments, removing mass less rapidly 
and thus overestimating the growth rate for super-Earths. For collisions among particles with 
these sizes, however, producing smaller fragments probably enhances collisional damping. If 
this enhanced collisional damping is sufficient to retain these small particles before low 
velocity destructive collisions remove them, super-Earth formation is likely. 

To summarize, fragmentation is a critical issue in the growth of super-Earths at 125- 
250 AU. For reasonable choices of parameters in Q* D , likely outcomes are similar to the results 
discussed in §4.1-4.2. Plausible choices for rrii have a broader range of consequences. In our 
assessment, the Morbidelli et al. (2009) algorithm probably prevents super-Earth formation; 
other algorithms (e.g., Kobayashi & Tanaka 2010; Weidenschilling 2010) probably allow 
super-Earth formation. Verifying these conclusions require additional simulations which are 
beyond the scope of the present effort. 


4-3.3. Initial Conditions 

Aside from the initial radius for planetesimals ro, our choices for the initial conditions 
have little impact on our results. Our adopted set of r Q spans a broad range of plausible 
sizes. Calculations with smaller r 0 will behave almost identically to simulations with r 0 = 
1 cm. Because collision times scale with initial radius, outcomes for swarms of planetesimals 
with ro > 1000 km will follow those for ro = 1000 km but on longer time scales. Thus, 
calculations with much smaller (larger) r 0 yield (do not yield) super-Earth mass planets. 

Another initial condition - the size distribution of planetesimals - is probably also not 
important. When planetesimals initially have a broad range of sizes, dynamical friction 
starts to modify the velocity distribution at the start of the calculation. Growth of large 
particles (r 0 = 1-100 km) then usually proceeds as much as a factor of two faster than 
growth from a mono-disperse set of particles (Kenyon & Bromley 2010, 2012). We expect 
similar differences at 125-250 AU. Thus, starting with a broad range of sizes might enable 
some initial conditions to yield super-Earths on 2-3 Gyr time scales instead of 4-6 Gyr time 
scales. 

A broad range of sizes probably has much less impact on the outcomes of calculations 
with small particles. When ro is small, dynamical friction and gravitational focusing are 
less important. Independent of the initial conditions, short collision times and slow growth 
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rapidly produce a broad range of sizes. Thus, our conclusions for the growth of small particles 
into super-Earths are fairly independent of the initial range of sizes. 

Our results indicate that the starting eccentricity also has little impact on the outcomes 
of the calculations. For simulations with r 0 > 1 km, our choices for the initial eccentricity 
span the likely range of possibilities. At smaller sizes, the equilibrium e and i are much 
smaller than our adopted eo and If interactions with the gas circularize the orbits of small 
particles effectively, the initial e and i might be much smaller. Gravitational instability of 
the swarm is then very likely (e.g., Michikoshi et al. 2007, 2009), leading to much more rapid 
growth of large particles. 


4-3.4- Dynamics 


The approximations used in coagulation codes begin to break down when (i) most of the 
mass is in a few planets and (ii) the dynamical interactions between these objects dominates 
collective interactions between small and large objects. To quantify these limits, we first 
consider dynamical constraints. At 125-250 AU, the escape velocity from the Solar System 
is roughly 2.5-4 km s -1 . With p s = 1.5 g cm -3 , planets have escape velocity 


R r 
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km s 
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Gravitational interactions raise the relative velocities of massless particles to roughly 60% 
of the escape velocity. Thus, planets with R p > 5000 km can eject small particles from the 
solar system. 

At the same time, collisional damping opposes gravitational scattering. Our calculations 
suggest damping overcomes scattering until planets reach super-Earth sizes. Thus, dynamical 
ejection of small particles is probably unimportant throughout most of each calculation. 

When protoplanets are large enough, their dynamical interactions are stronger and more 
chaotic than suggested by our Fokker-Planck approximation. Using the isolation criterion 
of eq. (2), two (four) equal-mass oligarchs on roughly circular orbits are isolated from one 
another when they have Rh ~ 0.028 a (0.014 a). These limits correspond to planets with 
masses of roughly 20 M® (3 M®). With typical masses of 2-3 M®, the two super-Earths 
produced in many of our calculations are clearly isolated. Large groups (10-20) of much 
lower mass planets are also stable. Although four super-Earths are nominally isolated, these 
planets are close to the dynamical stability limits and could be dynamically unstable (e.g., 
Kenyon & Bromley 2006; Bromley & Kenyon 2006, and references therein). 

Deriving the number of massive protoplanets in each calculations allows us to quantify 
the importance of strong dynamical interactions between protoplanets. Among 200 (150) 
calculations at 125 AU (250 AU), 45% (0%) produce one super-Earth, 47% (78%) produce 
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two, 6% (12%) produce three, and 2% (10%) produce four. In most of these calculations - 
98% at 125 AU and 90% at 250 AU - dynamical interactions among super-Earths are very 
unlikely. In others, dynamical interactions are probably common. 

During the other 800 (450) calculations at 125 AU (250 AU) which do not yield a 
super-Earth, protoplanets almost never interact dynamically. Because calculating these in¬ 
teractions is very cpu intensive (e.g., Kenyon & Bromley 2006), our neglect of these inter¬ 
actions is reasonable. Once we consider multiannulus calculations of super-Earth formation 
at 125-250 AU, we will revisit this issue. 

Although protoplanets eventually contain most of the mass remaining in each simulation, 
the approximations in our coagulation code rarely break down. Most simulations do not 
produce one or more super-Earths. In these situations, most of the mass is in protoplanets 
with sizes ranging from a few hundred km to a few thousand km. Our code accurately handles 
these situations. When 1-2 super-Earths form, these objects begin as 2000-3000 km objects 
which accrete material from the swarm of small particles with low eccentricity. Because 
these two objects do not interact dynamically, our statistical approach accurately derives 
the accretion rate of small particles. 

In a few cases, 3-4 oligarchs accrete small particles and grow to super-Earth masses. 
Our approach handles the growth of these oligarchs accurately. Dynamical interactions 
among these protoplanets become important near the end of the accretion process, when 
their masses approach an Earth mass. At this point our approximations begin to break 
down. Usually, we promote these objects into an n-body code and follow their dynamical 
interactions precisely (e.g., Kenyon & Bromley 2006; Bromley & Kenyon 2006, 2011a). Here, 
we simply note the formation of 3-4 objects with near super-Earth masses and save the more 
cpu intensive calculation for another study. 


5. DISCUSSION 

Our analysis indicates a plausible path for in situ production of super-Earth mass planets 
at 125-250 AU around solar-type stars. This path begins with a massive ring, Mo > 15 M 0 , 
composed of small particles, r 0 ~ 1 cm to 10 m. Successive mergers of these objects lead 
to larger and larger protoplanets. On a time scale of ~ 1 Gyr, coagulation yields several 
planets with masses m ~ 1-4 M 0 . 

Although we did not consider a range of initial masses for rings at 125-250 AU, previous 
studies demonstrate how the growth times and final masses of protoplanets depend on the 
initial mass and semimajor axis of a swarm of solid objects (e.g., Kenyon & Bromley 2008, 
2010). For all initial disk masses, the time scale for protoplanets to reach a fiducial radius 
(e.g., 1000 km) depends inversely on the initial disk mass. When Mo ~ 15 M 0 at 125 AU, 
it takes protoplanets 400 Myr to 1 Gyr to reach radii of 1000 km. Doubling (halving) the 
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initial disk mass halves (doubles) these time scales to 200-500 Myr (0.8-2 Gyr). Although 
the total masses of the most massive protoplanets also scale roughly linearly with mass, 
gravitational interactions among the protoplanets set the ultimate masses of the most massive 
protoplanets. For the conditions we consider, halving the initial disk mass reduces this 
ultimate mass by roughly 50%; however, doubling the initial mass may simply double the 
number of super-Earths rather than double the masses of the most massive super-Earths. 
Finally, the time scale for planet formation depends linearly on the orbital period. Doubling 
the initial semimajor axis nearly triples the formation times. 

Our conclusions raise several interesting issues concerning planet formation theory. 
Super-Earth formation at 100-250 AU also has direct observational consequences for the 
structure of debris disks around solar-type stars and for the architecture of the outer Solar 
System. After a brief discussion of several theoretical points, we describe how the new calcu¬ 
lations impact our understanding of debris disks and the outer Solar System. We conclude 
this section with several observational tests which can help to constrain plausible models. 


5.1. Theoretical Issues: Initial Conditions 

Our model for in situ super-Earth formation at 125-250 AU requires a massive reservoir 
of solids in the outer regions of the protosolar nebula. Observations of nearby star-forming 
regions suggest these initial conditions are plausible. In the p Oph and Taurus-Auriga molec¬ 
ular clouds, disks with radii of 100-300 AU are fairly common around pre-main sequence 
stars with ages of roughly 1 Myr (e.g., Andrews & Williams 2005; Andrews et ah 2009, 2010; 
Williams & Cieza 2011, and references therein). Of the twelve p Oph disks with outer radii 
larger than 100 AU, five have measured surface densities at least as large as the 0.02 g cm -2 
(125 AU) or the 0.005 g cm -2 (250 AU) required by our simulations. Within most disks 
around young stars, observations require particles at least as large as 1-10 mm (e.g., Williams 
& Cieza 2011). In many cases, observations also reveal that these particles have scale heights 
much smaller than the local scale height of the gas (e.g., Furlan et al. 2006, 2009). Thus, 
current evidence implies that these disks have particle sizes and total masses similar to those 
required by our calculations. 

Once cm-sized ‘pebbles’ form and settle towards the disk midplane, various processes 
shape their evolution. Collisions can slowly merge pebbles into larger objects or shatter them 
into fine dust grains (e.g., Birnstiel et ah 2010, 2012; Windmark et ah 2012; Garaud et ah 
2013). At 125-250 AU, viscous evolution slowly drives the gas (and any solids coupled to the 
gas) to larger distances (e.g., Lynden-Bell & Pringle 1974; Stepinski 1998; Chambers 2009). 
Large solids feel a headwind from the gas, which drags them towards (and perhaps into) the 
Sun (e.g., Adachi et ah 1976; Weidenschilling 1977; Rafikov 2004). However, the gas can also 
concentrate pebbles into much larger solid objects (e.g., Youdin & Goodman 2005; Johansen 
et ah 2007). Nearly all published studies of these phenomena focus on outcomes inside 



30 AU, where planet formation is a common process. Scaling these results to conditions at 
125-250 AU, it appears that 1 cm and larger particles can survive for the 1-3 Myr lifetime 
of the gaseous disk. Detailed studies of grain evolution are necessary to consider whether 
this conclusion is robust. If it is, our calculations suggest super-Earth formation is a likely 
outcome in large protostellar disks. 

Observations of solid objects beyond 30 AU might constrain the outer radius of the 
gaseous component of the protosolar nebula. Kretke et al. (2012) show that the gaseous disk 
excites Kozai oscillations in planetesimal leftovers from the formation of gas giant planets 
at 5-20 AU. Some of these planetesimals end up on orbits with low eccentricity and high 
inclination. The current lack of these objects suggests an outer edge of ~ 80 AU for the 
protosolar nebula. Observations over the next decade should allow significant refinement in 
this estimate. 


5.2. Theoretical Issues: Fragmentation 

The calculations in §4 illustrate the complex interplay between collisional damping, 
fragmentation, and the growth of the largest protoplanets in a gas-free environment. When 
protoplanets reach sizes of 500-1000 km, they initiate a collisional cascade which gradually 
grinds km-sized and smaller objects into micron-sized dust grains. Radiation pressure then 
ejects these grains from the planetary system. In our calculations, ejected grains comprise 
~ 25% to 80% of the initial mass in the ring. 

Throughout the collisional cascade, the largest protoplanets try to accrete smaller plan¬ 
etesimals and collision fragments. The rate these protoplanets grow depends on the magni¬ 
tude and timing of collisional damping (see also Kenyon & Bromley 2002a, 2004a; Goldreich 
et al. 2004). When the largest protoplanets and small fragments have relative eccentricity, 
€k,rei < 1 ( see ec f 9), protoplanets rapidly accrete the fragments. Super-Earth formation is 
likely. When either of these conditions is not met, protoplanets cannot accrete fragments as 
rapidly as the cascade grinds the fragments to dust. Protoplanets grow slowly and almost 
never reach super-Earth masses. 

In our calculations, the two conditions for super-Earth formation depend on the slope 
q of the power-law size distribution for collisional debris and the slope b in the relation for 
the mass of the debris produced in a collision ( m esc oc ( Qc/Q^b) b ] see eq. 6). At 125 AU, 
calculations with q & 4.5-5.5 and 6 ~ 9/8 are much more likely to meet these constraints 
on ek, re i than calculations with q & 3.5-4.0 and 6 « 1. At 250 AU, super-Earth formation is 
fairly independent of q and weakly dependent on b. 

Collisional damping is the physical mechanism that sets the final masses of protoplanets 
as a function of (q, b). When q is large, more fragments are concentrated into a smaller range 
of fragment masses. These fragments collide more rapidly, enabling collisional damping 
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to reduce collision velocities more effectively. Smaller velocities initiate a second runaway 
growth phase where protoplanets reach super-Earth masses. When q is small, fragments are 
spread out among a much larger range of masses. Damping is less effective; protoplanets 
grow more slowly. The b parameter magnifies this behavior. When b is large (small), high 
velocity collisions produce more (fewer) fragments, enhancing (reducing) collisional damping 
and the second runaway growth phase. 

Choosing the ‘best’ b and q for coagulation calculations is not straightforward. Despite 
recent advances in numerical simulations of debris production during binary collisions (e.g., 
Durda et al. 2004, 2007; Leinhardt & Stewart 2012; Benavidez et al. 2012), theory offers 
mixed guidance on several aspects of the algorithms used to model fragmentation. Currently, 
there is reasonable agreement on the binding energy Q* D in the gravity and strength regimes 
(see also Housen & Holsapple 1999, 2011). Modest changes consistent with the range of 
results from numerical simulations have modest impact on the growth of protoplanets (e.g., 
Kenyon & Bromley 2010, 2012). 

There is also reasonable agreement on the magnitude of b. For a broad variety of 
approaches, b = 1-1.25 (e.g., Davis et al. 1985; Housen & Holsapple 1990; Benz & Asphaug 
1999; Housen & Holsapple 1999; Kobayashi & Tanaka 2010; Leinhardt & Stewart 2012). For 
high velocity collisions between roughly equal mass objects, our approach is consistent with 
analytic derivations and numerical simulations. For cratering collisions where objects have 
very unequal masses, analytic fits to experimental data infer the ejected mass as a function 
of the impact velocity and the escape velocity of the more massive target (e.g., Housen & 
Holsapple 2011). When b = 1 or 9/8, our algorithm for the ejected mass generally agrees 
with these detailed fits (e.g., Kenyon et al. 2014). 

However, there is less agreement on r max the size of the largest fragment within the 
debris. For example, the recommended procedure of Leinhardt & Stewart (2012) results in 
two large fragments with total mass of roughly 10% of the total mass in the debris. From 
the results of the Durda et al. (2007) calculations, Morbidclli et al. (2009) derive a factor of 
5-10 smaller mass than Leinhardt & Stewart (2012). Our approach - where we place 20% 
of the total debris mass into the largest fragment - is similar to Leinhardt & Stewart (2012) 
but very different from Morbidelli et al. (2009). In future calculations, we plan to see if these 
differences have a clear impact on the overall outcomes of coagulation calculations. 

Moreover, there is little agreement on the size distribution for the debris. Leinhardt 
& Stewart (2012) propose a single power-law with q = 3.85 which extends from the largest 
remnant to the smallest particle allowed in the calculation. In contrast, Morbidelli et al. 
(2009) adopt q ~ 9-10 for the largest particles in the debris and q = 3.5 for the smallest 
particles. In both approaches, the largest particles contain most of the mass; thus, Morbidelli 
et al. (2009) have a much larger effective q than Leinhardt & Stewart (2012). Our calculations 
show that the magnitude of q has a dramatic impact on the size of the largest protoplanets. 
Using SPH/n-body simulations of binary collisions to derive a better appreciation of the 
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behavior of q as a function of the collision velocity, mass ratio, and composition would 
enable a better understanding of the outcomes of planet formation simulations. 

Finally, there has been little discussion of r min<d the minimum size of fragments produced 
during destructive collisions. Usually, r m i nd is zero; sometimes r min , td <C r maX)d . Together 
with the binding energy of planetesimals and larger protoplanets, r m i n _ d , r maXjd , and q es¬ 
tablish the importance of collisional damping and the likelihood of super-Earth formation 
in swarms of icy planetesimals. Robust predictions for the outcomes of planet formation 
calculations therefore require a broader discussion of these fragmentation parameters. 

In an interesting approach to these issues, Krijt & Kama (2014) infer a minimum size 
for debris fragments r min/ ] using constraints on the surface energy of colliding particles. In 
their model, a fraction / of the impact kinetic energy makes new surface area. When / = 
0, the surface energy of colliding particles is conserved. When / > 0, the collision produces 
additional surface area. For size distributions with 3 < q < 4, they derive r rmn , d as a function 
of the bulk properties of the colliding planetesimals. In high velocity collisions with fixed 
q , their picture produces size distributions similar to our algorithm. At lower velocities, 
collisions produce only a few (sometimes one) particles rather than a broad size distribution. 

Although the Krijt & Kama (2014) analytic model provides an interesting way to choose 
some aspects of the size distribution, it is valid only for 3 < q < 4 and does not directly con¬ 
strain / or q. Expanding on their methods might yield a more general analytic prescription 
for the size distribution of the debris and narrow the allowed range for q (see also Housen & 
Holsapple 1999). 

Without a useful analytic model, numerical collision calculations with a larger dynamic 
range in the sizes of debris particles would provide useful guidance on collision outcomes and 
the properties of the debris. Larger simulations would place better limits on plausible values 
for r mill)d and q in coagulation calculations. 


5.3. Applications: Debris Disks 

Aside from these theoretical issues, our results have clear consequences for the origin of 
the Solar System and debris disks at large distances from their host stars. Among known 
debris disks around solar-type stars, several have material at a ~ 100-200 AU. The disk 
in HD 107146 (G2 V spectral type, age ~ 100 Myr) extends from ~ 30 AU to ~ 150 AU 
(Corder et al. 2009; Ricci et al. 2014). Recent ALMA data indicate a reduction in dust 
surface density at 70-80 AU (Ricci et al. 2014). The large ring observed in HD 202628 
(G2 V, 2 Gyr) lies at ~ 200 AU and has a structure similar to the ring in Fomalhaut (Krist 
et al. 2012). In HD 207129 (G2 V, 2-3 Gyr; Krist et al. 2010; Marshall et al. 2011), a ring 
with roughly similar geometry is somewhat smaller (~ 160 AU). 
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Published analyses suggest fairly large masses of solids in HD 107146 and HD 207129. 
For HD 107146, fits to the dust distribution from ALMA data indicate a mass of roughly 
0.2 M® in cm-sized and smaller particles (Ricci et al. 2014). Deriving the total mass in solids 
requires an adopted size distribution. If n(r) oc r~ q for r — 1 cm to 100 km, likely values for 
q yield total masses ranging from 2-3 M® for q fa 4 to 100 M® for q ra 3.5 (see also Ricci 
et al. 2014). 

In HD 207129, detailed analysis of Herschel images and the spectral energy distribution 
implies a total mass of 0.007 M® in mm and smaller particles (Lohne et al. 2012). Adopting 
a power-law size distribution for 1 mm to 100 km objects yields masses ranging from 0.1 M® 
(q ~ 4) to 70 M® (q fa 3.5). 

Together with our numerical simulations, the large solid masses derived for these systems 
show that super-Earth mass planets are reasonable outcomes of the planet formation process. 
The geometries of the rings in HD 202628 and HD 207129 are similar to the geometry adopted 
for our rings at 125 AU and at 250 AU. With ages of 1-3 Gyr, our models suggest several 
super-Earths might be growing within the dusty debris of each system. 

In HD 107146, Ricci et al. (2014) infer a reduced surface density of dust at 70-80 AU 
within the extended disk of debris. If the lower surface density is a ‘gap’ in the disk produced 
by a planet, the likely mass of the planet is a few Earth masses. Based on our calculations, 
super-Earth formation at 70-80 AU around a 100 Myr old star requires a larger initial surface 
density than we considered for this study (see also Kenyon & Bromley 2008, 2010). With 
a factor of 3-5 increases in surface density, super-Earth formation in 100 Myr is a plausible 
outcome of our simulations (see also Ricci et al. 2014). 


5.4. Applications: Solar System 


The calculations in this study provide interesting constraints on the formation of a 
suspected planet X beyond 100 AU in the Solar System. Scaling the results at 125 AU and 
250 AU, the time scale to grow Earth-mass or larger planets from 1-10 cm pebbles is 
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where M 0 is the initial mass in an annulus of width 8a ~ 0.2a centered at semimajor axis a. 
Requiring planet formation to conclude within 3-4 Gyr allows super-Earths in annuli with 
Mo ~ 15 M® at a < 300 AU. More massive rings enable super-Earths outside 300 AU. 

Together with Bromley & Kenyon (2014), this study outlines plausible alternatives to 
standard migration models for super-Earths located at large distances from their host stars. 
In the migration hypothesis, a planet forms close to the host star and migrates out through 
a disk of gas or small planetesimals (e.g., Crida et al. 2009). In the scattering picture, a gas 
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giant planet dynamically scatters a much smaller planet from 5-10 AU to 100-200 AU; the 
gaseous disk then circularizes the orbit of the smaller planet at large a (Bromley & Kenyon 
2014). The in situ model envisions a massive swarm of solids in a ring at large a; coagulation 
within the ring allows the solid to grow into one or more super-Earth mass planets. 

In each of these scenarios, it is easier to understand a planet X in orbit at a ~ 100- 
250 AU than at a > 1000 AU. For the in situ picture, formation at a ~ 1000 AU requires a 
large reservoir of solids, > 300 M®. Fewer than 1% of the protostellar disks around Taurus- 
Auriga pre-main sequence stars and only 10% to 15% of the disks around Taurus-Auriga 
embedded protostars contain this much mass (see the discussion in Najita & Kenyon 2014, 
and references therein). In contrast, 20% to 60% of the disks around young stars contain 
enough mass to produce a super-Earth at 100-200 AU. Thus, in situ formation well beyond 
200-300 AU seems unlikely. 

Observations of young stars similarly limit the migration and scattering scenarios. In 
both mechanisms, the planet ends up on a roughly circular orbit within the gaseous disk 
surrounding the newly-formed central star. With typical characteristic radii of 20-200 AU 
(e.g., Andrews et al. 2010), disks around young stars rarely extend to 500-1000 AU. Thus, 
formation at 5-10 AU followed by migration/scattering to 500-1000 AU is also very unlikely. 


5.5. Tests 

Testing theories for super-Earths beyond 100 AU requires additional observations. Al¬ 
though recent observations identify Jovian mass planets associated with several circumstellar 
disks (e.g., Lagrange et al. 2010; Currie et al. 2014), direct imaging observations of isolated 
super-Earths with a > 100 ALI around other stars is unlikely with current or planned facili¬ 
ties (e.g., Hugot et al. 2012; Jovanovic et al. 2013; Macintosh et al. 2014). However, detecting 
super-Earths with surrounding debris clouds might allow direct tests as in Fomalhaut b (e.g., 
Kalas et al. 2005; Chiang et al. 2009; Kennedy & Wyatt 2011; Currie et al. 2012; Galicher 
et al. 2013; Kalas et al. 2013; Kenyon et al. 2014). LInless planet X has a moderate albedo, 
direct detection within the Solar System is also very challenging (Trujillo & Sheppard 2014). 
Observations of debris disks and Sedna-like objects provide alternate ways to test the various 
theoretical possibilities. 

Among debris disks, comparisons between observations and the debris predicted from 
collisional cascade or planet formation calculations allow indirect tests (e.g., Kenyon & Brom¬ 
ley 2002b, 2008, 2010; Lohne et al. 2012; Kobayashi & Lohne 2014). Although discrepancies 
between predictions and observations remain, these results generally confirm the cascade 
models (see also Gaspar et al. 2013). However, previous analyses did not include collisional 
damping of small particles within the cascade. In our calculations, damping enhances the 
population of 10 -2 — 10 2 cm particles and thus probably enhances the predicted emission 
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from small particles. As we learn more about the importance of damping as a function of 
the initial disk mass and semimajor axis, we will be able to compare the predicted emission 
with observations of complete ensembles of debris disks observed with IRAS, Spitzer, and 
Herschel. 

For nearby debris disks, ALMA observations offer the best hope of placing better limits 
on super-Earths at large distances from the parent star. When super-Earths form within a 
ring of small solid particles, they should clear out a gap with a half-width of roughly 3 Hill 
radii (e.g., Ida et al. 2000; Kirsh et al. 2009; Bromley & Kenyon 2011b, 2013). If ALMA 
detects narrow regions of low dust surface density in other disks and rings (as in HD 107146), 
comparisons with the predicted surface density within a gap might yield constraints on super- 
Earth formation in debris disks (e.g., Ricci et al. 2014). 

Within the Solar System, current evidence for a massive planet beyond 100 AU in the 
Solar System is limited and indirect. In the discovery paper for 2012 VP 113 , Trujillo & 
Sheppard (2014) suggest the need for a planet at 200-250 AU to explain the clustering of 
to among a dozen objects with perihelion distance q p > 30 AU and orbital eccentricity e > 
0.7. In addition to confirming the clustering, de la Fuente Marcos & de la Fuente Marcos 
(2014) note other clusterings in the orbital parameters. They propose two super-Earth mass 
planets at roughly 200 AU and 250 AU. However, constraints on the precession of the inner 
planets appear to constrain any planet X to much larger a > 1000 AU (e.g., Iorio 2009, 
2012, 2014, and references therein). Future studies need to reconcile the distribution of 
orbital parameters for Sedna-like objects with the precession rates for the gas giants and 
terrestrial planets. 

New calculations can place better constraints on the formation of planet X and lower 
mass objects like Sedna and 2012 VP 113 . For the in situ picture, ensembles of Earth-mass 
planets formed at 200 AU will eventually overcome collisional damping and scatter smaller 
objects (and possibly themselves) throughout the Solar System. Comparing the orbital 
elements of known objects on Sedna-like orbits with orbits derived from detailed numerical 
simulations of planet formation will allow good tests of this formation model (see also Kretke 
et al. 2012; Bromley & Kenyon 2014). 

Larger samples of Sedna-like objects can also test formation mechanisms. For 100- 
1000 km particles, coagulation calculations produce a size distribution with power-law slope 
q ph 4 (e.g., Kenyon & Bromley 2004d; Schlichting & Sari 2011; Kenyon & Bromley 2012; 
Schlichting et al. 2013). Among smaller objects, the slope has a broader range and depends 
on a variety of parameters. If 100-1000 km objects on Sedna-like orbits are as common as 
current surveys suggest (e.g., Trujillo & Sheppard 2014), measuring the slope of the size 
distribution could constrain the formation mechanism. 
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6. SUMMARY 

We analyze aspects of three mechanisms - migration, scattering, and in situ formation 
- for placing indigenous super-Earth mass planets on orbits with large semimajor axes, a > 
100 AU, around a solar-type star. In migration and scattering models, a super-Earth forms 
within a few Myr at a ~ 5-10 AU. Interactions with other protoplanets and a massive 
disk then stabilize the orbit at large a. The in situ picture requires a massive reservoir of 
solid material where super-Earths form directly. For all three hypotheses, super-Earths at 
100-300 AU are much more likely than super-Earths beyond 300 AU. 

Using a suite of single-annulus coagulation calculations, we investigate the physical 
requirements to form super-Earth mass planets at 125-250 AU around solar-type stars. 
In annuli with Mq & 15 M©, super-Earth formation is common. When this mass is in 
planetesimals with ro ~ 1 cm to 1 m, super-Earths form in 1-3 Gyr at 125 AU and in 
2-5 Gyr at 250 AU. Although swarms of larger planetesimals (r 0 ~ 10 m to 10 km) can 
produce super-Earths in < 4 Gyr at 125 AU, growth at 250 AU takes longer than the age of 
the Solar System. Thus, in situ growth of super-Earths at 200-300 AU requires a massive 
ring of small planetesimals. 

Collisional damping of the debris from destructive collisions is an essential ingredient in 
super-Earth formation at 125-250 AU. During the first 300 Myr to 1 Gyr of evolution, the 
usual pattern of slow growth, runaway growth, and oligarchic growth leads to the formation 
of 1000-3000 km protoplanets and a swarm of leftover small (0.1-10 km) planetesimals 
with large relative velocities (see also Kenyon 2002; Kenyon & Bromley 2002b, 2004a, 2008, 
2010). Destructive collisions among these leftovers generates a multitude of nun-sized to 
m-sized fragments. Damping reduces the relative velocities of the fragments by 1-2 orders 
of magnitude, to ek, re i < 1 where ek, re i = e/r^ and th = (mmax/SMo) 1 / 3 is the relative 
size of the Hill sphere of the most massive protoplanet. The small fragment velocities enable 
a second runaway growth phase, where protoplanets grow rapidly. When protoplanets have 
£k,rei < U they reach super-Earth sizes in 1-2 Gyr. 

Fragmentation also plays a crucial role in the formation of super-Earths at 125-250 AU. 
In our calculations, the debris from destructive collisions has a power-law size distribution, 
n{r) oc r~ q . At 125 AU, setting q = 4.5-5.5 promotes the formation of super-Earths. When q 
= 3.5-4.0, super-Earth formation is very rare. At 250 AU, all q m 3.5-5.5 allow super-Earth 
formation. In models with q fa 4.5-5.5, however, super-Earths tend to be more massive 
than in calculations with q fa 3.5-4.0. Within the fragmentation algorithm, the variation 
of ejected mass with impact energy is also important. In our approach, m esc oc ( Q c /Qf)) b ■ 
Calculations with b = 9/8 tend to produce more massive protoplanets than calculations with 
b = 1. 

Improving the calculations requires a better understanding of several parameters in the 
fragmentation algorithm. In our approach, the mass of the debris ejected during a collision 
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is m esc = 0.5 (Q c /Qr)) b , where Q c is the collision energy, Q* D is the collision energy required 
to eject half of the mass involved in the collision, and b = 1.00-1.25. For the mass of the 
debris, we set uil — 0.2 m esc and adopt a constant value q for slope of the power-law size 
distribution of the debris. In a future study, we plan to investigate the impact of adopting 
a functional form for the largest object in the debris, m,L = rriL,o(Qc/Qb)~ bL where ~ 
0.5-1.0 (e.g., Kobayashi et al. 2010; Weidensc-hilling 2010). However, analytic studies and 
more comprehensive numerical simulations of binary collisions would enable more informed 
choices for 0 , b, bL, and q as a function of the collision energy and the bulk properties of 
the planetesimals. 

Observations of debris disks allow interesting tests of the models. Expanding our single- 
annulus models to a full multi-annulus calculation will enable comparisons between predic¬ 
tions and observations of the time-dependent debris disk luminosities. ALMA observations of 
the dust surface density in the gaps of debris disks can place limits on super-Earth formation 
rates. 

Discovery of planet X and additional dwarf planets on Sedna-like orbits would also test 
the plausibility of in situ super-Earth formation at 125-250 AU. From the mass distribution 
of protostellar disks, indigenous super-Earths with a > 500 AU are unlikely (see also Najita 
& Kenyon 2014). Robust estimates of the semimajor axis of any planet X directly constrain 
planet formation theory. Aside from placing better constraints on an orbit for planet X 
(Trujillo & Sheppard 2014; dc la Fuente Marcos & de la Fuente Marcos 2014), orbital elements 
and (eventually) the size distribution for new members of the inner Oort cloud can test the 
predictions of the coagulation models. 


We acknowledge generous allotments of computer time on the NASA ‘discover’ cluster. 
We thank the referee, K. Ohtsuki, for a constructive and thorough review. Advice and 
comments from T. Currie and M. Gcllcr also greatly improved our presentation. Portions of 
this project were supported by the NASA Astrophysics Theory and Origins of Solar Systems 
programs through grant NNX10AF35G and the NASA Outer Planets Program through grant 
NNX11AM37G. 


A. Appendix 

Orchestra contains many algorithms to follow the formation and evolution of a plan¬ 
etary system. With few analytic solutions available, verifying these algorithms requires 
comparisons with other approaches. In previous papers (e.g., Kenyon & Luu 1998; Kenyon 
& Bromley 2001; Bromley & Kenyon 2006; Kenyon & Bromley 2008; Bromley & Kenyon 
2011 a), we compare our numerical results with analytic solutions and n-body simulations. 
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Here, we describe several additional tests. 


A.l. Runaway and Oligarchic Growth 


Accurate calculations for runaway and oligarchic growth depend on the spacing factor 
6 = mk+i/rrik between adjacent mass bins (e.g., Ohtsuki & Nakagawa 1988; Ohtsuki et al. 
1990; Wetherill & Stewart 1989; Wetherill 1990). When bins have fixed masses, adopting 
large 6 risks severe underestimates in the growth times for massive objects. Ohtsuki et al. 
(1990) derive the constraints on 5 which allow accurate solutions (see al;so Fraser 2009). In 
our approach, mass bins have fixed boundaries but variable masses (see also Wetherill & 
Stewart 1993; Kenyon & Luu 1998). Calculations with 5 > 1 always lag a ‘perfect’ solution, 
reaching a fiducial size later in time. Larger <5’s always produce larger lags. 


To understand the origin of the lag, we consider a simplified example. With preset 
boundaries rrik+ 1/2, the average mass of objects in a bin is fhk, where rrik+1/2 < hik+i < iTik+ 3/2 
and rrik+3/2 = 5 • mk+i/2- When an object in bin k + n merges with an object in bin k, the 
mass of the new object is m c = fhk + rhk+n■ In an ideal calculation with an infinite number of 
mass bins, the mass of the new object is placed in a new bin which is distinct from all other 
bins. In a real calculation, however, placing this object in a new bin with index k c > k + n 
requires m c > m,k+ n + 1/2- To derive a simple estimate of when this inequality is satisfied, we 
set fhk+n = S n fhk■ The constraint on promoting the new object to a new bin then becomes: 


Q S n + 1 )m k 

fin+l/2 mk 


(Al) 


Throughout a calculation, fhk is distributed uniformly in the range (j -1 / 2 -^ 1 / 2 m k ■ The 
minimum fhk which satisfies eq. (Al) is 5~ l / 2 m k - We can then write a much simpler inequality 
for promotion to a new bin: 

S n+1 - S n - 1 < 0 . (A2) 

When this inequality is satisfied, promotion to a new bin always occurs. Otherwise, promo¬ 
tion fails some fraction of the time. This failure is the source of the lag in all coagulation 
calculations with variable mass bins. The fraction of failures grows with increasing 5. Thus, 
the lag grows with larger 5. 

Before describing results of our test calculations, it is worth examining several exam¬ 
ples. When n = 1, we consider a simple collision between objects in adjacent mass bins. 
Calculations with 5 < 1.62 satisfy the inequality in eq. (A2). In solutions with 5 > 1.62, 
collisions fail to promote the merged object into a higher mass bin roughly 20% of the time. 
When merged objects are not promoted, the average mass of the largest object in the grid 
is smaller than the mass of a promoted object. Smaller objects have smaller collision rates 
and longer collision times. Thus, calculations with <5 > 1.62 always lag an ideal calculation 
for this simple merger. 
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As n increases, we consider collisions between objects with larger and larger mass ratios. 
Automatic promotion to the next bin requires smaller and smaller 5. When n = 3, 10, 25, or 
62, promotions never fail when 5 < 1.38, 1.19, 1.10, 1.05. For collisions with large n, smaller 
6 yields smaller failure rates. Thus, the lag declines with decreasing 5. 

Previous discussions of the magnitude of the lag as a function of 5 focus on the evolution 
of the overall size distribution (Wetherill 1990; Kenyon & Luu 1998; Morbidelli et al. 2009). 
Here, we illustrate the lag by following the growth of the largest object. In the first example, 
we consider growth when the cross-section is Aij = 77/7777 where 7 is a constant. Defining 
the dimensionless time 77 = 'jriot where no is the initial number of particles with mass mo, 
there is an analytic solution for the size of the largest object as a function of no and rj 
(Wetherill 1990). 

Fig. 17 compares results for four values of 5 with the analytic solution. Independent of 
5, r max follows the same evolution: slow growth for rj ~ 0.0-0.95 and exponential growth for 
77 a; 1. The precise timing of exponential growth depends on 5. Models with 5 = 2 reach 
the runaway when rj = 1.091, lagging the exact solution by 9.1%. When 5 = 1.05, runaway 
begins when rj = 1.001. With a lag of only 0.1%, these calculations almost precisely match 
the analytic solution. 

Although solutions for A l} oc rn^rrij provide excellent tests of a coagulation code, realistic 
cross-sections for particles in circumstellar disks are not as extreme (see also Wetherill 1990; 
Lee 2000; Malyshkin & Goodman 2001). When the collision rate depends solely on the 
geometric cross-section, A^ oc (+ mj) 2 ^ 3 . If gravitational focusing is important, A l} oc 
(m, + rrij) 4 / 3 . Calculating the evolution for either cross-section yields a more appropriate 
assessment of the lag as a function of 5 in actual protostellar disks. With no analytic solution 
available, we adopt calculations with 5 = 1.10 as the reference. For clarity, we rescale the 
dimensionless time to place exponential growth for 5 = 1.10 at rj = 1. 

Fig. 18 compares results for six values of 5 with A VJ oc (777 + m-k ) 4 / 3 (see also Wetherill 
1990). All solutions have the same form. Compared to the analytic model in Fig. 17, the 
largest object grows somewhat faster for 77 = 0.0-0.95; close to 77 = 1, runaway growth is 
somewhat less explosive. The timing of the runaway is very sensitive to 5. Models with 5 = 
2 reach runaway roughly 40% later in time than models with 5 = 1.05. 

When A^ oc (m* + nrij ) 2 / 3 , growth is insensitive to 5 (Fig. 19). Solutions with <5 = 2 lag 
calculations with 5 = 1.1 by less than 2 %. Calculations with 5 = 1.25-1.7 have lags ranging 
from 0.5-1.5%. Thus, mass spacing has little impact on growth with purely geometric cross- 
sections (Wetherill 1990). 

These examples suggest that 5 is an important factor in models of runaway growth 
but not slow or oligarchic growth (see also Wetherill 1990). To confirm this point, Fig. 20 
compares the evolution of r max for complete evolutionary calculations with different values 
of 5. During an early phase of slow growth (when t < 1 Myr), 1 cm particles reach maximum 
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sizes of a few meters on time scales that are independent of 8. After runaway growth begins, 
small particles rapidly grow to km sizes. Once particles have r max « 0.1 km, calculations with 
larger 5 enter the oligarchic phase before those with smaller 5. Eventually all calculations 
reach oligarchic growth; the rate of growth is then fairly independent of 5. After several Gyr 
of evolution, all models have r max « 3000-7000 km (Fig 20). 

When 0.1 km < r rnax < 100 km, the growth times in these calculations depend on 5. 
For 5 < 1.2, the time scale to produce 100 km objects is t wo & 500-800 Myr. Although 
calculations for each 5 exhibit a range in f 10 o, the median time to reach r max ~ 100 km 
is independent of 8. For larger 8, the median growth time lags the 8 < 1.2 solutions by 
10%-20% (<5 = 1.41) to 25%-50% (<5 = 2.0). Thus, calculations with 8 < 1.2 yield a better 
estimate of the time scale to produce 100 km objects. 

As objects grow larger than 100 km, the growth time depends less on 8 and more on 
the fragmentation history. For any 8, collisional damping of small fragments can enable a 
second runaway growth phase at t ~ 100 Myr to 1 Gyr. The timing and magnitude of this 
phase depends on (i) the amount of material available in small objects and (ii) the ability 
of collisional damping to reduce particle velocities significantly. Sometimes conditions allow 
the production of massive objects with r max > 10 4 km. Other simulations produce objects 
with r max fa 3000-4000 km. 

A.2. Collisional Damping and Gravitational Stirring 

In our coagulation calculations, we derive changes in particle velocities from collisional 
damping, dynamical friction, gas drag, and viscous stirring. Kenyon & Luu (1998, Appendix) 
and Kenyon & Bromley (2008, §2) describe the algorithms for each of these processes. Col¬ 
lisional damping is a central issue in this study; thus, we summarize our treatment for a 
single annulus. We then discuss results for several standard examples. 

To treat collisional damping in a single annulus, we solve two sets of differential equa¬ 
tions 2 : 

fjh 2 l=k 

—-jp = ~y mih * ~ mkh ^ ~ mk + Wki) (A3) 

1=0 

and 

—77^ = "ir ("W - rn k v\ - (m k + mi)vl) Ii(P k i) (A4) 

at 1=0 Pkl 


2 This discussion is adapted from the multiannulus approach of eqs. 8-9 from Kenyon & Bromley (2008). 
For a single annulus, we drop the subscripts for annulus i and annulus j and eliminate the term which treats 
the overlap between adjacent annuli. 
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where C in = a cM e ki pi V kt F gM (■ r k + n) 2 , fih = (: i\ + tf)/(el + ef), and pi is the volume 
density of planetesimals with mass mi (Wetherill & Stewart 1993; Kenyon & Luu 1998). In 
these expressions, a co u = 0.57-0.855 (Wetherill & Stewart 1993), e jk = 1/2 for k = l and 1 
for k ^ l, V k i is the relative particle velocity, and F gkl is the gravitational focusing factor. 
The integrals I e and are elliptic integrals described in previous publications (e.g., Hornung 
et al. 1985; Wetherill & Stewart 1993). 

Together with similar equations for coagulation, dynamical friction, gas drag, and vis¬ 
cous stirring (see Kenyon & Luu 1998; Kenyon & Bromley 2002a, 2008, and references 
therein), we solve these equations explicitly using an adaptive time step which conserves 
mass and kinetic energy to machine accuracy. We use more accurate and more cpu intensive 
fourth-order Runge Kutta and Richardson extrapolation algorithms to verify the explicit 
solutions. Kenyon & Luu (1998) summarize comparisons of our single annulus solutions for 
terrestrial planet formation with those of Wetherill & Stewart (1993). 

To demonstrate our ability to treat collisional damping accurately, we consider several 
idealized problems. In the first test, we follow Levison & Morbidclli (2007) and consider 
the evolution of 1 cm particles in a single annulus at 30-35 AU. The particles have initial 
eccentricity eo = 0.1, inclination i 0 = eo/2, mass density p s = 1 g cm -3 , and total mass 
Mo = 10 M®. The calculations include gravitational stirring and collisional damping, but 
particles cannot merge or fragment during a collision. 

Our results closely match those summarized in Levison & Morbidelli (2007) and Mor- 
bidelli et al. (2009). Particle random velocities damp rapidly and reach an equilibrium at 
e eq = 9.3 x 10 -10 and i eq = 5.1 x 10“ 10 in 7300 yr (Fig. 21). The e-folding times of t e ~ 
350 yr and t t ~ 400 yr are close to the Levison & Morbidelli (2007) results of 320 yr and 
360 yr. The ratio tjti = 0.875 is essentially identical to their t e /ti = 0.889. 

To compare this approach to analyses of Saturn’s rings (e.g., Goldreich & Tremaine 
1978; Bridges et al. 1984; Supulver et al. 1995; Porco et al. 2008), we modify the damping 
algorithm. Defining the coefficient of restitution as 

cr = c 0 v;^ R , (A5) 

with c 0 = 0-1 and kn = 0.1-0.3, the kinetic energy loss rate is 5E ig = A^v^l — cj^/Vij where 
Vij is the relative velocity of two particles and V tJ is the volume occupied by the particles. 
Replacing our standard loss term with SE tJ allows us to compare damping rates derived with 
two different approaches. 

The grey lines in Fig. 21 plot results for c 0 = 0.99 (upper curve) and c 0 = 0.01 (lower 
curve) with k ki = 0.15. hi these solutions, the e-folding times for the eccentricity are t e ~ 
435 yr (c 0 = 0.99) and 425 yr (c 0 = 0.01); for inclination, t t m 435 yr (c 0 = 0.99) and 425 yr 
(co = 0.01). Standard collisional damping in our calculations proceeds on time scales similar 
to models with a very small coefficient of restitution cq < 10~ 3 (see also Kenyon & Luu 1999, 



40 


and references therein). 

To test the evolution in an annulus relevant to the calculations described in the main 
text, we examine damping and stirring in a ring of mono-disperse particles with p s = 
1.5 g cm -3 and Mq = 15.79 M® at 112.5-137.5 AU. For the damping calculations, we 
set eo = 10 2 e eq and io = eo/2. To address time scales for gravitational stirring, we set 
e 0 = 10~ 2 e eq and i 0 = eo/2. 

In both sets of simulations, particles reach an equilibrium where the 3D velocity dis¬ 
persion is roughly 60% of the escape velocity. The equilibrium has h/v = 0.49, i/e = 0.55, 
and e eq = 2.38 x 10~ 4 (f?/l km). The time scales to reach equilibrium are typically much 
longer at 125 AU than at 30-35 AU (Figs. 22-23). When gravitational stirring dominates 
(Fig. 22), particles reach equilibrium on time scales ranging from ~ 1000 yr for 1 cm particles 
to roughly 1 Gyr for 10 km particles. Damping time scales are roughly a factor of 30 longer 
(Fig. 23). 


Aside from validating the algorithms for stirring and damping, these simulations help 
us to set likely initial conditions for e and i at 125-250 AU. For a swarm of particles with 
velocity dispersion c p , surface density E, and angular velocity D, the swarm is gravitationally 
unstable when c p VL < 7tGE (e.g., Chiang & Youdin 2010). Setting E = E 0 (a/1 AU) _3//2 , the 
swarm is gravitationally stable when the eccentricity has e > e s , where 


10 


-4 


30 g cm 


-2 


\ 1/2 


100 ALU 


(A6) 


For our initial surface density at 125-250 AU, the minimum e for gravitational stability is 
e s « 10" 4 . 


Particles with r > 0.4 km have e eq larger than e s . At 125-250 AU, these particles 
reach equilibrium on time scales longer than the lifetime of the protosolar nebula and other 
protoplanetary disks (e.g., Dauphas & Pourmand 2011; Williams & Cieza 2011; Cloutier 
et al. 2014, and references therein). Thus, we can safely set eo somewhere between e s and 

^■eq- 

Although particles with radii r ~ 10-100 m have e q < e s , the time scale to reach 
this equilibrium is longer than the disk lifetime. On these time scales, particles also grow 
considerably (see Figs. 1-2). Thus, we set e 0 ~ e a . 

Particles with much smaller radii r ~ 1-100 cm reach very small equilibrium e on 
relatively short time scales, ~ 0.05-5 Myr. For our adopted E and eo = 10 -4 , these particles 
grow on similar time scales. To avoid the complications of gravitational instabilities, we also 
set e 0 = e s for these particles. 

In any swarm of particles at 125-250 AU, the likely e and i depend on the evolutionary 
history, which is sensitive to the structure of the disk. In this paper, the adopted eo and % 
are intermediate between the e eq and values expected for a turbulent disk. The calculations 



in §4 describe how these choices impact our results. The discussion in §5 considers how to 
improve these choices. 
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Fig. 1.— Growth of the largest object at 125 AU as a function of initial planetesimal size 
for calculations with b = 9/8 and q = 3.5. Although simulations starting with smaller 
planetesimals produce large planets faster than those with larger planetesimals, all systems 
with Mo ~ 15.8 M 0 and r o < 10 km yield planets with r max m 3000-6000 km. In rings with 
r 0 > 100 km, collisions are too infrequent to produce planets with r max > 500-2000 km. 
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Fig. 2. 


As in Fig. 1 for calculations with 6 = 1. 
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Fig. 3.— As in Fig. 1 for calculations with r 0 = 1 cm or 100 m, b = 9/8, and various q as 
indicated in the legend. In calculations with q > 4.5, planets reach super-Earth masses in 
1-3 Gyr. 
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Fig. 4.— Maximum planet radius at a = 125 AU as a function of ro and q for b = 9/8. 
To improve clarity, points are slightly offset from the nominal r 0 . The horizontal grey bar 
indicates the radius of an Earth-mass planet with p s = 1.5 g cm -3 . Calculations with r 0 = 
1 cm to 1 km and q > 4.5 produce super-Earth mass planets. 
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Fig. 5 


As in Fig. 4 for calculations with 6 = 1. 
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Fig. 6.— Comparison of cumulative (lower panel) and relative cumulative (upper panel) size 
distributions for three idealized power law size distributions. To construct the relative size 
distribution, we adopt a normalization with q s = 2.7 which yields n c , re i = 1 at r = 1 cm. 
Black curves: rti — 1, q s — 2.7, qi = 1.0, qi = 3.5, r\ = 0.1 km, r 2 = 10 km. Cyan curves: ni 
— 1, q s — 2.8, qi = 1.2, qi = 4.0, r\ = 0.03 km, r 2 = 30 km. Magenta curves: ni — 1, q s — 
2.6, qi = 1.7, qi = 4.0. r\ = 0.3 km, r 2 = 100 km. Relative size distributions allow an easier 
assessment of the differences between power-law components. 
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Fig. 7.— Evolution of the relative cumulative size distribution n c _ re i for a calculation with 
r 0 = 1 cm, B = 9/8, and q = 3.5. Here, n C)Te i = n c /n r where n c is the cumulative number 
distribution from the calculation, n r = nor -2,5 , and n o = 10 23 . Numbers to the left of each 
curve indicate time in Myr. 
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Fig. 8. 


As in Fig. 7 for later epochs. 
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Fig. 9.— Evolution of the relative eccentricity eu in Hill units for a calculation with ro = 
1 cm, B = 9/8, and q = 3.5 at times (in Myr) indicated in the legend. As the maximum 
particle size grows, small particles maintain e# ~ 10. At larger sizes (r > 0.03 r max ), the 
eccentricity follows a power law eu ~ 0.2(r/r max ) 1 ■ 
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Fig. 10.— As in Fig. 9 for later epochs. As the evolution proceeds, collisional damping 
depresses the eccentricity for small particles. Once these particles have been swept up or 
destroyed, the eccentricity of small particles rises. 
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Fig. 11.— Relative cumulative size distributions when r max = 3000 km for calculations with 
q = 3.5 (black points), 4.5 (green points), and 5.5 (cyan points). Calculations with larger q 
have more material in particles with r — 1 cm to 0.1 km and less material in particles with 
r < 1 cm. 
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Fig. 12.— As in Fig. 11 for the relative eccentricity. In calculations with larger q, particles 
with r > 1 cm have smaller ek t rei] smaller particles have larger e^^ei- 
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Fig. 13.— Growth of the largest object as a function of initial eccentricity. For the two 
particle sizes considered (ro = 10 cm and 1 km), the legend in the lower right corner indicates 
the initial eccentricity relative to the nominal eo- Overall, growth is fairly independent of 
the initial e and i. 
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Fig. 14.— Growth of the largest object at 250 AU as a function of q for b = 1 and r 0 
= 1 cm and 10 m. For rings of small particles with r 0 = 1-10 cm, coagulation produces 
super-Earth mass planets for all q in 2-4 Gyr. Larger particles with ro > 10 m cannot grow 
into super-Earths over the age of the Solar System. 
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Fig. 15.— Maximum planet radius at a = 250 AU as a function of r 0 and q for b = 1. 
To improve clarity, points are slightly offset from the nominal r^. The horizontal grey bar 
indicates the radius of an Earth-mass planet with p s = 1.5 g cm -3 . Calculations with ro 
= 1 cm to 1 m produce super-Earth mass planets. Larger planetesimals evolve into much 
smaller planets on 1-10 Gyr time scales. 
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Fig. 16.— As in Fig. 15 for calculations with b = 9/8. On 1-4 Gyr time scales, super- 
Earths form from planetesimals with r 0 = 1 cm to 1 m. Occasionally, ensembles of 10 m 
planetesimals can grow into 1-2 super-Earths if q > 4.5. Larger planetesimals produce 
planets with radii ranging from a few hundred km to a few thousand km. 
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Dimensionless Time 


Fig. 17.— Mass of the largest object as a function of dimensionless time rj for cross-section 
Aij oc miirij. The solid black line plots the analytic solution. Colored points show results 
for numerical solutions with different values for the mass spacing factor 5 as summarized in 
the legend. The numerical solutions lag the analytic solution by 9.1% (5 = 2.0), 5.7% (5 
= 1.7), 2.5% (5 = 1.4), 2.0% (5 = 1.25), 0.3% (5 = 1.1), and 0.1% (5 = 1.05). Although 
calculations with 6 = 1.05 follow the analytic result almost exactly, solutions with 5 < 1.4 
yield satisfactory approximations to the exact result. 
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Fig. 18.— As in Fig. 17 for Ay oc (rn t rrij )A' 3 . The numerical solutions lag the 5 = 1.10 
result by 43.7% (5 = 2.0), 37.4% (5 = 1.7), 31.1% (5 = 1.4), 22.5% (5 = 1.25), and 9.0% (5 
= 1.17). 
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Fig. 19.— As in Fig. 18 for A,j oc ( rriimj ) 1 ^ 3 . The numerical solutions lag the 5 = 1.10 
result by 1.9% (5 = 2.0), 1.6% (5 = 1.7), 1.3% (5 = 1.4), 1.0% (5 = 1.25), and 0.5% (5 = 
1.17). 



68 



Time (yr) 


Fig. 20.— Growth of 1 cm particles at a = 125 AU as a function of mass spacing parameter 
<5. Although all calculations produce objects with similar r max , simulations with large <5 lag 
those with smaller 6. 
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Fig. 21.— Collisional damping within a ring of small particles at 30-35 AU (Levison & 
Morbidelli 2007). The ring contains 10 M e in 1 cm particles. Heavy lines: damping for 
inelastic collisions. Light lines: damping for elastic collisions with coefficient of restitution 
Cq = 0.99 (upper line) and Cq = 0.01 (lower line). 
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Fig. 22.— Gravitational stirring within a ring of particles at 112.5-137.5 AU. The ring 
has surface density £ = 0.0215 g cm -2 and contains 15.8 M e in mono-disperse particles. 
Solid and dashed lines: evolution of the eccentricity e relative to the equilibrium eccentricity 
eo = 2.38 x 10 _4 (i?/l km) for 1 cm, 10 cm, 1 m, 10 m, 100 m, 1 km, and 10 km particles. The 
time to reach e _1 of the equilibrium eccentricity is t s ti r = 1.67(72/1 km) Myr. Vertical grey 
bar: epoch of planetesimal formation and gas dispersal in a protoplanetary disk (e.g., Klcine 
et al. 2009; Dauphas & Pourmand 2011; Williams & Cieza 2011, and references therein). 
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Fig. 23.— As in Fig. 22 for collisional damping. The e-folding time is t ( [ amp = 
10 9 (i?/l km) yr. 







